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ABSTRACT 


During  high-speed  injection  in  small-diameter  nozzles,  a  combustible  mix¬ 
ture  may  experience  high  enough  temperature  for  premature  ignition  to  oc¬ 
cur.  Among  the  various  physical  mechanisms  which  may  lead  to  such  an 
undesirable  effect,  shear-induced  heating  is  believed  to  play  an  important 
role.  The  objective  of  this  effort  is  to  study,  theoretically  and  numerically, 
the  impact  of  shear-induced  heating  on  the  likelihood  of  mixture  ignition. 

The  manisfestation  of  shear-induced  heating  effects  is  believed  to  be  con¬ 
siderably  dependent  on  the  mode  of  injection.  For  the  purpose  of  this  study, 
two  injection  modes  (short-duration  injection  and  continuous  injection)  are 
distinguished.  For  short  duration  injection,  shear-induced  heating  is  expected 
to  be  confined  to  thin  laminar  boundary  layers.  This  situation  is  analyzed 
using  computational  codes  which  implement  simplified  physical  models  for 
high-speed  injection  in  a  constant-diameter  nozzle.  The  computational  codes, 
which  numerically  integrate  the  incompressible  vorticity  transport  and  energy 
equations,  are  applied  to  predict  the  maximum  temperatures  experienced  by 


11 


the  mixture.  The  impact  of  injection  speed  and  inlet  mixture  temperature, 
and  the  effect  of  variable  mixture  properties  were  also  investigated.  For  con¬ 
tinuous  injection,  boundary  layer  transition  is  expected  to  occur.  In  order 
to  characterize  the  shear  heating  mechanisms  within  a  transitional  or  turbu¬ 
lent  flow  environment,  one  must  characterize  the  mean  flow  and  temperature 
field,  and  fluctuations  around  the  mean.  To  this  end,  the  mean  tempera¬ 
ture  solutions  are  constructed  using  well-established  empirical  correlations 
for  the  mean  velocity  field.  Meanwhile,  the  impact  of  temperature  fluctu¬ 
ations  is  evaluated  in  light  of  direct  simulation  of  isotropic  turbulence  and 
transitional  channel  flow. 
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Chapter  1 
Introduction 


One  of  the  concerns  in  design  and  operation  of  nozzles  used  for  the  injec¬ 
tion  of  combustible  mixtures  is  prevention  or  at  least  minimization  of  hazards 
associated  with  preignition  of  the  mixture  inside  the  nozzle.  Such  an  undesir¬ 
able  event  may  damage  or  even  cause  a  complete  failure  of  the  entire  system. 
In  particular,  the  behavior  of  high-speed  liquid  propellant  flow  has  a  direct 
bearing  on  safety  issues  (Knapton  et  al.  1992).  Among  the  various  mech¬ 
anisms  that  could  lead  to  this  undesirable  event,  it  is  believed  that  viscous 
heating  and  compression  of  entrained  gas  bubble  (Nigmatulin  k  Khabeev 
1974;  Bourne  k  Field  1991;  Field  1992;  Yuan  and  Prosperetti  1996)  play  an 
important  role.  Our  attention  is  focused  here  on  shear-induced  (or  viscous) 
heating  effects. 

Shear-induced  heating  effects  arise  due  to  the  rapid  acceleration  of  the 
mixture  (or  fluid)  to  sufficiently-high  velocities  which  are  required  to  deliver 
the  necessary  charge  to  a  combustion  chamber  (or  mixing  device).  Thus, 
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high  shearing  rates  are  established  which  may  result  in  high-enough  temper¬ 
atures  for  spontaneous  combustion  to  occur. 

The  manifestation  of  shear  induced  heating  effects  are  believed  to  be  con¬ 
siderably  dependent  on  the  mode  of  injection.  For  the  purpose  of  this  study, 
we  distinguish  between  the  following  two  injection  modes:  (i)  short-duration 
injection,  in  which  the  charge  is  impulsively  accelerated  through  the  nozzle, 
and  (ii)  continuous  injection,  where  a  steady  mass  flow  rate  develops  within 
the  nozzle. 

To  illustrate  the  differences  between  the  two  modes,  we  consider  the  fol¬ 
lowing  typical  range  of  nozzle  parameters  and  operating  conditions.  For  both 
cases,  the  characteristic  nozzle  diameter  and  length  are  taken  as  2  mm  and 
20  mm,  respectively,  while  the  characteristic  mixture  velocity  is  roughly  200 
m/s-  The  mixture  properties  are  roughly  the  same  as  those  of  water,  so 
that  its  density  and  kinematic  viscosity  are  respectively  approximated  by 
p  =  lOOOfcy/m^  and  u  =  10“®m^/s.  For  short-duration  injection,  the  pulse 
duration  is  assumed  to  be  approximately  10ms. 

For  continuous  injection,  the  flow  field  dynamics  -  and  corresponding 
shear  heating  rates  -  may  be  appropriately  interpreted  in  terms  of  the  flow 
Reynolds  number.  In  the  targeted  range  of  nozzle  applications  estimated 
above,  the  Reynolds  number  based  on  the  nozzle  diameter  is  in  the  order 


of  10®.  Accordingly,  both  transitional  and  laminar  boundary  layers  are  ex¬ 
pected  to  develop.  In  the  transitional  case,  shear-induced  heating  is  a  result 
of  “laminar”  heating  within  the  boundary  layer,  and  intermittent  rollup  of 
“large-scale”  vortices.  When  the  boundary  layer  becomes  turbulent,  shear- 
induced  heating  is  determined  by  the  mean  turbulent  temperature  distribu¬ 
tion  (including  the  “laminar”  heating  within  the  viscous  sublayer),  and  the 
fluctuating  component  (which  is  dominated  by  the  presence  of  wide  spectrum 
of  eddy  structures). 

In  the  first  part  of  this  study,  it  is  argued  that  the  characterization  of  the 
flow  field  as  a  fully-developed  channel  flow  (whether  laminar,  transitional, 
or  turbulent)  is  not  appropriate  for  short-duration  injection.  In  this  case, 
both  the  evolution  flowfield  and  shear-induced  heating  mechanism  must  be 
interpreted  in  terms  of  the  thin,  unsteady,  developing  boundary  layers  which 
are  generated  along  the  nozzle  walls.  To  verify  this  claim,  we  start  by  esti¬ 
mating  the  typical  thickness  of  the  boundary  layers  due  to  the  diffusion  of 
the  vortex  sheets  which  form  due  to  the  impulsive  acceleration  of  the  mix¬ 
ture.  The  scaling  of  the  diffusion  zone  (or  boundary  layer)  may  be  conserva¬ 
tively  estimated  through  the  following  measure  of  the  displacement  thickness 
d*  (i/t)^/^  =  0.1mm.  Thus,  during  the  injection  period,  viscous  effects  re¬ 
main  confined  to  a  small  region  close  to  solid  boundaries.  Most  of  the  flow  in 
the  channel  is  still  unaffected  by  viscous  diffusion  from  the  wall  and  remain 
potential.  Noting  that  the  maximum  Reynolds  number  based  on  displace- 


ment  thickness  is  less  than  2  x  10^,  and  considering  the  stabilizing  effects  of 
the  acceleration  of  the  potential  core,  it  is  highly  unlikely  that  any  instability 
mechanism  resulting  in  a  transition  to  turbulence  amplifies  significantly  dur¬ 
ing  such  injection  durations.  Consequently,  shear-induced  heating  associated 
with  pulsed  injection  modes  must  be  analyzed  using  the  unsteady,  primarily 
laminar,  development  of  thin  boundary  layers  close  to  the  nozzle  walls. 

Based  on  the  above  scaling  arguments,  it  is  obvious  that  different  ap¬ 
proaches  are  necessary  while  simulating  shear-induced  heating.  For  continu¬ 
ous  injection,  the  risk  of  mixture  ignition  may  be  best  analyzed  by  investi¬ 
gating  the  following  two  possible  scenarios: 

•  The  mean  temperature  distribution  exceeds  the  ignition  threshold  in  a 
region  of  the  flow. 

•  The  mean  temperature  distribution  is  everywhere  below  the  ignition 
threshold,  but  sustained  temperature  excursions  of  large-enough  am¬ 
plitude  may  locally  exceed  this  “critical”  value. 

In  analyzing  the  latter  possibility,  we  are  less  concerned  with  high-frequency, 
high-amplitude  turbulent  fluctuations,  since  these  fluctuations  are  least  effec¬ 
tive  in  sustaining  high  shear  rates,  and  are  not  considered  as  a  likely  precur¬ 
sor  to  mixture  ignition.  Thus,  the  contribution  of  these  fluctuations  may  be 
modeled  as  part  of  the  mean  “turbulent”  temperature  distribution.  However, 
it  is  conceivable  that  (intermittent)  rollup  of  large-scale  energetic  eddies  in  a 


transitional  boundary  layer  may  sustain  large-enough  viscous  heating  rates 
for  the  temperature  to  exceed  to  the  critical  value. 

However,  as  noted  before,  it  is  believed  that  the  above-mentioned  phe¬ 
nomena  are  highly  unlikely  during  short-duration  injection.  In  this  case,  the 
primary  concern  is  the  maximum  temperature  attained  at  or  near  the  nozzle 
boundaries.  Thus,  the  physical  and  numerical  modeling  efforts  must  focus 
on  obtaining  accurate  estimates  of  the  wall  and  boundary-layer  temperature 
distributions.  The  likelihood  of  non-Newtonian  behavior  should  also  be  ex¬ 
plored  for  this  type  of  flow. 

Following  the  above  arguments,  we  shall  first  focus  on  a  short-duration 
acceleration  of  an  incompressible  Newtonian  fluid  in  a  constant-diameter  noz¬ 
zle.  Specifically,  we  were  concerned  with  a  simplified  channel  geometry  which 
still  enables  us  to  observe  all  of  the  essential  shear  heating  mechanisms,  and 
aimed  at  determination  of  peak  temperatures  experienced  by  the  liquid  un¬ 
der  flow  conditions  where  the  nozzle  wall  boundary  layer  is  predominantly 
laminar. 

In  order  to  obtain  conservative  estimates  of  peak  temperatures,  an  un¬ 
steady  one-dimensional  analysis  of  the  temperature  distribution  associated 
with  sudden  fluid  acceleration  over  a  flat  insulated  boundary  is  first  con¬ 
ducted  (Chapter  3).  The  analysis  yields  analytical  expressions  of  the  pre- 


vailing  temperature  distribution  during  the  early  stages  of  boundary  layer 
formation.  In  particular,  the  obtained  expressions  relate  the  peak  wall  tem¬ 
perature  to  the  prevailing  Eckert  and  Prandtl  numbers.  Results  reveal  a 
quadratic  dependence  of  the  normalized  wall  temperature  on  impulse  ve¬ 
locity,  and  a  complex  nonlinear  variation  with  Prandtl  number.  The  latter 
dependence  highlights  a  concern  for  mixtures  having  high  Prandtl  number. 

Next,  simulation  of  high-speed  flow  in  an  axisymmetric  nozzle  are  per¬ 
formed  (  Chapter  4).  The  numerical  schemes  are  based  on  a  finite-diflference 
discretization  of  vorticity-based  momentum,  and  energy  equations.  Variants 
of  the  formulation  are  considered  which  use  one  of  the  following  approaches: 

1.  full,  unsteady  equations  of  motion  on  a  rectangular  grid 

2.  unsteady  parabolized  equations  of  motion  on  uniform  mesh 

3.  unsteady  parabolized  equations  on  stretched  grid 

4.  steady  parabolized  equations  on  stretched  grid 

A  numerical  study  of  all  four  formulations  is  performed  and  used  to: 

•  establish  the  validity  of  adopting  a  parabolized  approximation 

•  validate  computational  results  and  reflne  grid  stretching  techniques 

•  numerically  analyze  the  impact  of  inlet  velocity  and  temperature  bound¬ 
ary  conditions. 


Finally,  the  effects  of  wall  heat  transfer  (Chapter  5)  and  temperature- 
dependent  properties  (Chapter  6)  are  discussed. 

In  the  second  part  of  this  study,  we  are  focused  on  modeling  of  shear- 
induced  heating  mechanisms  for  continuous  injection  periods,  where  a  steady 
mass  flow  rate  develops  within  the  channel.  For  continuous  injection,  flow- 
field  dynamics  in  the  entire  channel  section  are  important.  Accordingly, 
one  must  be  able  to  characterize  the  transitional  or  turbulent  flow  dynamics 
across  the  channel,  i.e  the  mean  flow  and  fluctuations  about  the  mean  flow. 

In  Chapter  7,  we  obtain  the  estimates  for  the  mean  flow  using  simple 
empirical  correlations  for  the  mean  velocity  and  energy  dissipation,  which 
yields  an  analytical  expression  of  mean  temperature  solution.  Next,  we  have 
tried  to  characterize  the  temperature  fluctuations  of  forced  isotropic  tur¬ 
bulent  flow  using  direct  numerical  simulation  (in  Chapter  8).  At  last,  the 
direct  numerical  simulation  is  performed  on  turbulent  channel  flow  (Chapter 
9),  where  transitional  turbulent  boundary  layer  is  observed,  and  the  impact 
of  turbulent  fluctuations  is  examined. 


Part  I 

Numerical  Study  of  Laminar 

Flow 
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Chapter  2 

Numerical  approach 


2.1  General 

As  mentioned  in  the  introduction,  we  are  primarily  interested  in  charac¬ 
terizing  shear-induced  effects  during  high-speed  injection  of  a  combustible 
mixture  into  a  small  diameter  nozzle. 


For  such  applications,  shear-induced  heating  effects  are  primarily  con¬ 
centrated  in  the  viscous  boundary  layers  which  develop  close  to  the  nozzle 
boundaries.  Furthermore,  since  these  effects  depend  strongly  on  shear  rates, 
they  will  be  primarily  manifested  in  the  thin  section  of  the  nozzle,  where  the 
mixture  is  accelerated  to  high  velocities.  Therefore,  a  proper  resolution  of 
boundary  layer  development  in  the  thin  section  of  the  nozzle  is  clearly  seen 
to  be  a  crucial  step  in  any  analysis  aiming  at  the  characterization  of  peak 
temperatures  experienced  by  the  mixture. 
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In  general,  accurate  resolution  of  the  temperature  within  the  thin  bound¬ 
ary  layers  developed  in  the  thin  nozzle  section  requires  a  complex  analysis 
which,  in  particular,  must  account  for  the  flowfield  dynamics  upstream  of 
the  nozzle  section.  This  is  the  case  because  the  flow  conditions  at  the  inlet 
of  the  nozzle  may  be  strongly  dependent  on  the  geometry  of  the  converging 
section  and  on  the  upstream  flow  profile. 

The  potential  differences  in  inlet  conditions  from  different  nozzle  geome¬ 
tries  may  be  illustrated  by  considering  a  smoothly  tapered  nozzle,  and  a 
nozzle  formed  by  a  sudden  contraction.  No  flowfield  separation  is  expected 
in  the  first  case,  resulting  in  a  smooth  and  monotonic  boundary  growth  across 
the  contraction.  When  the  contraction  is  sudden,  on  the  other  hand,  flowfield 
separation  may  occur  upstream.  In  this  situation,  the  boundary  layer  within 
the  thin  nozzle  section  may  develop  starting  from  a  minute  or  vanishingly 
small  thickness.  Since  the  largest  shear  rates  scale  according  to  the  ratio  of 
the  characteristic  velocity  in  the  potential  core  (U)  to  the  inlet  boundary 
layer  thickness  (5o)j  shear-induced  heating  mechanisms  are  expected  to  be 
significantly  more  pronounced  in  the  latter  case. 

The  potential  dependence  of  shear  heating  rates  on  nozzle  geometry  and 
upstream  flow  conditions  complicates  numerical  study  of  the  correspond¬ 
ing  heating  effects.  On  one  hand,  the  simulation  must  tackle  the  geometric 
complexity  of  the  contraction  section  -  whether  smooth  or  sudden  -  and 


accurately  describe  inflow  conditions.  On  the  other  hand,  detailed  flowfleld 
resolution  will  necessitate  restricting  the  analysis  to  a  limited  class  of  nozzle 
geometries  and  flow  conditions. 

In  order  to  overcome  these  limitations,  a  simplified  approach  to  the  prob¬ 
lem  will  be  adopted  in  the  first  part  of  this  study.  The  approach,  which  is 
motivated  by  our  desire  to  minimize  the  likelihood  of  mixture  ignition,  calls 
for  conservatively  estimating  peak  temperatures  by  focusing  on  worst-case 
scenarios.  To  this  end,  the  analysis  shall  assume  that  the  flow  conditions  at 
the  inlet  of  the  thin  nozzle  section  correspond  to  a  vanishingly  small  boundary 
layer  thickness  and  an  essentially  flat  velocity  profile.  By  doing  so,  simula¬ 
tion  of  flow  dynamics  upstream  of  the  nozzle  is  avoided,  and  upper  bounds 
on  the  mixture  temperature  can  be  efficiently  obtained  as  a  function  of  the 
inlet  velocity,  the  nozzle  diameter,  and  the  mixture  properties. 


2.2  Formulation 

The  physical  formulation  used  in  the  present  study  is  based  on  the  following 
assumptions: 

•  The  combustible  mixture  is  an  incompressible  Newtonian  liquid. 


•  The  flowfleld  within  the  nozzle  remains  axisymmetric. 
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Until  section  6,  we  shall  furthermore  assume  that  the  mixture  has  con¬ 
stant  properties,  so  that  its  motion  is  governed  by  the  continuity,  momentum 
and  energy  conservation  equations,  respectively  written  as: 


V-  u=  0 

(2.1) 

Dv  1 

^ =  vp  +  V 

L)t  p 

(2.2) 

DT 

Dt  ^  + 

(2.3) 

where  v  is  the  velocity  vector,  p  is  density,  p  is  pressure,  1/  is  the  kinematic 
viscosity,  Cp  is  specific  heat  at  constant  pressure,  k  is  the  thermal  con  ductiv- 
ity,  $  is  the  viscous  dissipation  function,  and  ^  v  -A  is  the  material 

derivative.  For  an  axisymmetric  flowfield,  the  governing  equations  are  recast 
as: 


dvr  Vr  dv, 

+  (2.4) 

dVr  dVr  dVr  1  dp  fd'^Vr  I  dVr  Vr  d^Vr.  ,  , 

+  +  +  +  (2.5) 


and  the  viscous  dissipation  function  is  given  by: 

T  rfcr/^^r\2  /  \2  fdVr  dv^^n 

®  =  (^.s) 

In  the  numerical  study  of  high-speed  nozzle  flow,  we  shall  also  rely  on 
the  vorticity  form  of  the  momentum  equations.  This  alternative  formulation 
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is  obtained  by  taking  the  curl  of  the  momentum  equations  and  using  the 
continuity  equation,  to  get: 


du 

dt 


+  + 


d  ,  .  d'^u  d  .u). 


(2.9) 


where 


(2.10) 


is  the  vorticity.  The  latter  is  related  to  the  radial  and  streamwise  velocity 
components  through  a  scalar  streamfunction,  V",  such  that: 


1  dx}) 

r  dz 
1  dip 

r  dr 


Substituting  the  velocity  components  from  Eq.  (2.10)  into  Eq.  (2.9)  results 
in  the  familiar  vorticity  streamfunction  relationship: 
d  dip.  1  d'^'^ 

dr^r  dr^  ^  r  dz'^  ^  (2-11) 

The  advantage  of  using  a  vorticity-streamfunctionn  formulation  in  axisym- 
metric  flow  is  that  the  continuity  equation  is  naturally  satisfled  since  the 
velocity  fleld  is  expressed  as  the  curl  of  a  vector  potential  (ip/r),  and  the 
number  of  “unknown”  fleld  variables  is  reduced  from  three  in  the  primitive 
variables  formulation  {vr,  v^,  p)  to  two  in  the  present  case  {ip,  u). 


2.3  Normalization 


The  numerical  simulations  discussed  in  the  following  sections  solve  a  normal¬ 
ized  form  of  the  governing  equations.  This  normalized  form  is  obtained  by 
nondimensionalizing  variables  with  respect  to  the  appropriate  combination 
of  the  mixture  density,  p,  the  nozzle  radius,  i2,  the  maximum  velocity  at  the 
nozzle  inlet,  U,  and  the  mixture  inlet  temperature,  To.  For  the  choice  of 
characteristic  density,  length,  velocity  and  temperature  scales,  the  normal¬ 
ized  constant-property  governing  equations  are  listed  in  Table  2.1  below. 
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Table  2.1 

Normalized  Model  Equations  for  a  Constant  property 
Mixture  in  Axisymmetric  Coordinates 


Governing  Equations 
Vorticity  transport 

dt  +  +  a;*?)  + 

Streamfunction 


_a  1^  lay 

dr  r  dr  ^  r  dz^  ^ 


Energy 


dT  ^  dT  dT  1  Ec 

dt  dr 


Definitions 

Velocity  field 


Vr  =  -^^ 

V  = 
z  r  dr 


Viscous  dissipation 


$  =  2  ,  /'^\2  ,  ('^\2]  ,('^r  dVz  5 


Parameters 
Reynolds  number 


(2.12) 


(2.13) 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


Prandtl  number 


(2.18) 


Eckert  number 


Chapter  3 

Impulsively-Started  Quasi-ld 
Flow 


The  model  equations  given  in  the  previous  section  indicate  that,  in  the 
present  formulation,  the  prevailing  flow  regime  within  the  nozzle  is  described 
in  terms  of  three  governing  parameters:  the  Reynolds  number  Re,  the  Prandtl 
number  Pr,  and  the  Eckert  number  Ec.  In  order  to  gain  an  appreciation  for 
the  impact  of  these  parameters  on  shear  induced  heating  mechanisms,  we 
first  examine  the  simplified  problem  of  impulsively  started  motion  over  a  flat 
insulated  plate. 

This  physical  setting,  which  closely  approximates  the  early  stage  of  bound¬ 
ary  layer  formations  at  the  nozzle  walls,  is  specified  as  follows.  The  flat  plate 
has  infinite  extent  and  coincides  with  the  a:-axis  in  the  2D  Cartisian  x  -  y 
plane.  The  fluid  is  at  rest  for  times  t  <  0,  and  is  impulsively  accelerated  to  a 
velocity  U  at  t  =  0.  The  flow  is  unbounded  in  the  +y  direction  and  remains 
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one-dimensional.  Under  these  conditions,  the  dimensional  constant-property 
momentum  equations  reduce  to: 
du  B^u 

It  ~ 

with  boundary  conditions: 

fu(!/  =  0,i)  =  0 
I  u{y  -4  -foo,t)  =  U 

Here,  u  denotes  the  streamwise  velocity  component,  and  the  initial  con¬ 
dition  may  be  simply  expressed  as:  u{y,t  =  0)  =  0. 

The  solution  of  the  above  heat  equation  for  the  streamwise  velocity  com¬ 
ponent  is  given  by  the  well-known  similarity  solution[l]: 


u{r])  =  Uerf{r}) 

(3.3) 

where  rj  =  y/ y/Aut  is  the  similarity  variable,  and 

erf{ri)  =  -7=  /  exp{-x^)dx 

Jo 

(3.4) 

denotes  the  error  function. 

Meanwhile,  the  dimensional  form  of  the  energy  equation  takes  the  form: 

dT  d^T  V  ,du,o 

dt  ^ dy"^  Cp  dy^ 

(3.5) 

with  boundary  and  initial  conditions  respectively  given  by: 

[  f(?/  =  o,t)  =  o 

\  T{y  -^00)  =  To 

(3.6) 

II 

II 

(3.7) 
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3.1  Solutions  of  the  Energy  Equation 

Analytical  solutions  of  the  energy  equation  (3.5-3. 7)  are  sought  using  two 
different  techniques.  First,  we  rely  a  Green’s  function  approach  which  is 
based  on: 


(i)  enlarging  the  solution  domain  to  the  entire  plane  — oo  <y  <  -Foo 

(ii)  defining  the  following  extension  to  the  velocity  field: 

»•(!;,()  =  t/er/(i)  (3,8) 

(iii)  exploiting  the  symmetry  of  the  extended  velocity  field  by  recasting  the 
energy  equation  as: 


dT*  d^T*  u  .du\o  d^T*  . 

di  dy^  Cpdy  dy^  -KCpt  2ut 

with  boundary  and  initial  conditions  respectively  given  by: 


(3.9) 


T*{y  ±00,  t)  =  To 
T*{y,t  =  0)=To 


(3.10) 

(3.11) 


An  analytical  expression  for  the  time-dependent  solution  to  the  above 
equations  is  readily  expressed  using  the  Green's  function  of  the  heat  equation; 
we  have: 


X  exp( 


TTCp  Jo  Y  T(2Q;(t  —  r)  +  ut) 
,2 

)dT 


y 


2{2a{t  —  r)  +  ut)' 


(3.12) 


While  evaluation  the  above  integral  is  generally  cumbersome,  simple  expres¬ 
sions  for  the  wall  temperature,  =  T{y  =  0)  =  T*{y  =  0)  may  be  readily 
obtained.  The  analytical  results  are  summarized  by: 


7-  +  ¥  Pr  =  2  (3.13) 

I  1  +  f  \/E H]/Pr(Pr  -  2)  +  Pr  -  1)  Pr  >  2 

It  is  interesting  to  note  that  the  wall  temperature  is  time-independent 
and  the  above  expressions  do  not  explicitly  depend  on  the  Reynolds  num- 
ber,due  to  the  absence  of  a  characteristic  flow  length  scale. 


Pr  <2 
Pr  =  2 


In  order  to  obtain  solutions  for  the  temperature  profile,  and  to  check  on 
the  validity  of  the  above  expressions,  a  similarity  variable  approach  is  also 
implemented.  To  this  end,  we  introduce  a  new  similarity  variable. 


P 


=  y 

y/Ui 


(3.14) 


and  postulate  that  the  temperature  distribution  is  a  function  of  P  only,  i.e. 
T  =  fiP)-  Substituting  this  into  Eq.  (3.5)  yields: 


(3.15) 


where  primes  refer  to  diflferentiation  with  respect  to  p.  Multiplying  Eq. 
(3.15)  by  the  integrating  factor  jAa)  and  collecting  terms  gives: 


Integrating  Eq.  (3.16)  once  with  respect  to  ^  and  enforcing  the  adiabatic 
wall  constraint  /'(O)  =  0  yields: 

exp(£/?^)/'  =  [--!>(- -  £))d.  (3.17) 

In  general,  Eq.  (3.17)  is  integrated  numerically.  However,  for  special  val¬ 
ues  of  Prandtl  number,  the  integration  may  be  easily  performed  analytically. 
Specifically,  for  Pr  =  1,  we  have: 

m  =  (3.18) 

and  in  normalized  form: 

=  1  +  ^(1  -  erf{r]f)  (3.19) 

Meanwhile,  for  Pr  =  2,  we  get: 

21  P  (P 

f{(5)=To  +  —  exp(-^)  (3.20) 

TT  Cp  Zi 

and  in  normalized  form: 

T{r})  2  ^  /  9x 

-7j^  =  ^  +  -Ecex.^{-7f)  (3.21) 

It  is  easily  verified  that  the  wall  temperature  predicted  using  the  above  ex¬ 
pressions  coincides  with  that  obtained  using  the  Green’s  function  approach. 


3.2  Discussion 


Results  of  the  above  theoretical  analysis  are  depicted  in  Figs.  (3. 1-3.3).  Fig. 
(3.1)  shows  normalized  temperature  profiles  for  mixtures  with  different  Eck¬ 
ert  and  Prandtl  numbers.  In  Figs.  (3.2)  and  (3.3),  a  mixture  having  the 
same  heat  capacity  as  water  and  initially  at  29b°K  is  assumed;  normalized 
and  dimensional  temperature  distributions  are  plotted  for  different  injection 
velocities  and  different  Prandtl  numbers. 

Figs.  3.1  —  3.3  enable  us  to  clearly  visualize  the  roles  of  the  Eckert  and 
Prandtl  numbers  on  shear  heating  during  the  early  stages  of  boundary  layer 
formation.  In  particular,  Fig.  3.1  clearly  reflects  the  linear  variation  of  the 
normalized  wall  temperature  with  Eckert  number,  while  Figs  3.2  -  3.3  reflect 
a  quadratic  dependence  of  the  impulse  velocity.  In  addition,  the  nonlinear 
dependence  of  peak  temperature  on  the  Prandtl  number  is  also  illustrated. 
The  logarithmic  divergence  of  the  wall  temperature  as  Pr  ->  oo  underscores 
a  particular  concern  for  mixtures  having  high  Prandtl  number. 

Finally,  it  is  interesting  to  note  that  the  shear  heating  predictions  of  the 
present  ID  analysis  are  only  dependent  on  the  Eckert  and  Prandtl  num¬ 
bers.  Thus,  an  explicit  dependence  on  the  fluid  density  is  lacking  from  the 
present  predictions.  Accordingly,  the  manifestation  of  shear-induced  heating 
mechanisms  is  expected  to  be  similar  for  incompressible  flow  conditions  char¬ 
acterized  by  the  same  Eckert  and  Prandtl  numbers,  even  if  the  corresponding 


mixtures  have  significantly  different  densities. 


Ec=0.001 


Ec=0.004 


Ec=0.01 


Ec=0.04 


Figure  3.1:  Normalized  temperature  profiles,  plotted  against  the  similarity 
coordinate  77,  for  different  Eckert  numbers  and  Prandtl  numbers.  The  value 
of  the  Eckert  number  is  indicated. 


Figure  3.3:  Heating  of  a  mixture  impulsively  accelerated  over  a  flat  plate. 
The  mixture  is  assumed  to  have  the  same  heat  capacity  as  water  and  an 
initial  temperature  of  2%°K.  Temperature  profiles  are  plotted  against  the 
similarity  coordinate  77,  for  different  injection  speeds  and  Prandtl  numbers. 
The  value  of  injection  velocity  is  indicated. 


Chapter  4 

Axisymmetric  Flow 


We  now  turn  our  attention  to  the  numerical  modeling  of  high-speed  injection 
in  an  axisymmetric  nozzle.  Throughout  this  section,  we  shall  assume: 

(i)  an  incompressible  mixture  with  constant  properties, 

(ii)  that  the  nozzle  has  radius  =  2mm  and  length  L  =  20mm,  and 

(iii)  that  the  nozzle  walls  are  insulated. 

Thus,  following  the  discussion  of  Section  2,  we  are  interested  in  the  simulation 
of  the  normalized  vorticity  transport  and  energy  equation,  respectively: 


m  +  +  3:(r)  + 


Re^dr^  '  dr^r'  '  dz^ 

dT  ^  dT  ^  dT  1 

—  +  Vr^  +  v^—  =  — — v^r  +  — $ 
or  oz  RePr  Re 


dt 


(4.1) 

(4.2) 


The  boundary  conditions  associated  with  the  above  system  of  governing 
equations  are  summarized  as  follows: 
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1.  At  the  nozzle  inlet,  the  velocity  and  temperature  distributions  are  as¬ 
sumed  to  be  known.  Furthermore,  we  assume  that  the  radial  velocity 
component  vanishes  at  inflow,  so  that  the  velocity  boundary  conditions 
is  given  in  terms  of  a  streamwise  velocity  profile,  Ui{r).  The  latter  is 
used  to  derive  streamfunction  and  vorticity  boundary  conditions  re¬ 
spectively  through: 

V'inW  =  ^  C«z(C)c?C  (4.3) 

and 


y  V  /  \ 

(4.4) 

Thus,  inflow  boundary  conditions  are  expressed  as: 
r  ^(r,2  =  0)  =  '0m(r) 

<  cj(r,2:  =  0)  =a;i„(r)  (4.5) 

[  T{r,z  =  0)  =  Tin{r) 


2.  At  the  nozzle  wall,  the  potential  and  no-slip  velocity  conditions  are  im¬ 
posed,  as  well  as  the  adiabatic  surface  boundary  condition.  We  thus 
have: 


i’{r  =  l,z)=f^ru^{r)dr 
|^(r  =  l,z)  =  §:(r  =  l,2)  =  0 


(4.6) 


3.  At  the  centerline,  axial  symmetry  conditions  are  enforced,  and  the  stream- 
function  value  is  arbitrarily  set  to  zero.  Accordingly,  we  have: 

f  ^(r  =  0,  z)  =  0 
I  f{r  =  0,z)  =  0 


(4.7) 


4.  At  the  downstream  computational  boundary,  derivative  outflow  condi¬ 
tions  are  imposed;  we  have 

|^(r,  2  =  20)  =  ^(r,  z  =  20)  =  ^(r,  z  =  20)  =  0  (4.8) 

4.1  Inflow  and  Initial  Conditions 

Since  the  fine  details  of  the  mixture  injection  into  the  nozzle  are  generally  not 
known  or  difficult  to  determine,  impulsively-accelerated  flow  should  be  ide¬ 
ally  considered.  Unfortunately,  unlike  the  ID  theoretical  study,  this  situation 
is  not  easily  accommodated  in  the  context  of  a  multi-dimensional  numerical 
simulation.  First,  the  early  stages  of  growth  of  viscous  boundary  layers, 
which  initially  form  as  infinitely-thin  vortex  sheets,  cannot  be  adequately 
captured  using  finite  grid  sizes.  Second,  the  imposition  of  a  flat  inlet  veloc¬ 
ity  profile  would  unnecessarily  burden  the  computations,  since  these  cannot 
handle  the  associated  singular  behavior  near  the  nozzle  wall. 

In  order  to  overcome  these  difficulties,  the  initial  velocity  filed  is  assumed 
to  admit  a  small,  but  finite-thickness  boundary  layer  at  the  nozzle  wall.  The 
velocity  profile  corresponds  to  the  analytical  solution  given  in  the  previous 
section,  evaluated  at  a  time  significantly  smaller  than  the  injection  duration. 
Moreover,  the  imposed  inflow  velocity  profile  is  also  taken  to  coincide  with 
the  same  boundary  layer  solution.  A  “brute  force”  approach  is  then  adopted 
which  calls  for  decreasing  the  initial  boundary  layer  until  the  computed  so¬ 
lution  is  essentially  independent  of  this  “free  parameter” . 
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Meanwhile,  two  different  approaches  are  followed  in  specifying  the  initial 
temperature  distribution  and  the  inlet  temperature  boundary  condition.  The 
first  is  analogous  to  that  used  for  the  velocity  profile  and  is  based  on  using 
the  ID  solution  to  appropriately  fit  a  thermal  boundary  layer  at  the  nozzle 
wall.  The  second  approach  simply  assumes  a  flat  inlet  temperature  profile. 
As  shown  in  the  computations  discussed  below,  both  approaches  yield  nearly 
identical  peak  temperature  predictions. 


4.2  Direct  Numerical  Simulation 


Direct  numerical  simulation  of  the  governing  equations  is  first  performed. 
This  approach  is  based  on  accurate  and  detailed  resolution  of  all  relevant 
length  and  time  scales  using  the  full  equations  of  motion.  To  this  end,  the 
vorticity  transport  and  energy  equations  are  first  discretized  in  time  using 
the  3rd-order  Adams  Bashforth  scheme,  while  the  viscous  dissipation  term 
is  treated  using  the  2nd-order  Crank-Nicolson  approach.  Thus,  the  time- 
discretized  vorticity  and  energy  equations  are  expressed  as: 


(4.9) 

(4.10) 


where 
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_  d(uvr)  d(uv^) 

~  dr  ^  dz 

(4.11) 

„  _  d'^u)  d  Lo  d^u 

dr^  r  dz^ 

(4.12) 

^  dT  dT 

R  =  Vr-^  +  V^  — 
dr  dz 

(4.13) 

(4.14) 

respectively  denote  the  convective  vorticity,  vorticity  diffusion  convective 
temperature  and  thermal  diffusion  terms,  while  (70, 71, 72)  =  (23/12,  -16/12, 5/12) 
are  the  integration  constants  of  the  3rd-order  Adams-Bashforth  scheme. 

Spatial  discretization  of  the  above  equations  is  performed  on  regular  rect¬ 
angular  grid  of  mesh  size  (Ar,  A2:).  All  “internal”  spatial  derivatives  are 
treated  using  second-order  centered  differences,  resulting  in  the  expressions 
summarized  in  Table  4.1.  Second-order  treatment  of  temperature  boundary 
conditions  at  the  wall  is  performed,  but  one-sided  first-order  differences  are 
used  at  outflow  and  the  centerline.  The  standard  first-order  treatment  of  the 
vorticity  boundary  condition  at  the  wall  is  also  used. 
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Table  4.1 

Discretized  Form  of  Derivative  Terms  in  Standard  Notation 


Vorticity  Convection 

2Ar 

2Az 


Pi 


V 


In 

+- 


Vorticity  Diffusion 

n..  cv  -  ^^i,j  +  ^t-1.7  ,  COi+ij-Ui-ij  Wij 

^^3  ^  A  I  ^ 


Ar^ 
A«2 


2r,Ar 


Temperature  Convection 


i?jj  iVr)i,j  2Ar 

Thermal  Diffusion 

S  ! 


2Az 


Tj+ij  2Ttj  +  Tj-ij  ^  1  Pj+ij  ~  Tj-ij 


Ar2  '  n  2Ar 

PjJ+l  ~  ^^,3  +  Pi,j-1 


(4.15) 


(4.16) 


(4.17) 


Az^ 

Viscous  Dissipation 

$  off  ~  I  (Mi,js^2  ,  (i'^z)i,j+l  -  {Vz)i,j- 

2Ar  '  r.-  ^  oA-y 

-1  -  h 

2Az  '  2Ar 


(4.18) 


2Ar  '  •  '  r,  '  ■  '  2Az 

I  /  (^r)ij+l  —  {Vr)i,j-1  (^2)i+l,j  ~  (l’z)i-lj-.2 

^  2Az  2Ar  ^ 


1\2 


(4.19) 
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Vorticity-Streamfunction 

'‘Pi+lj  ~  2'0t,j  +  _  1  V’i+l.j  ~  fpi-lj 

Ar^  n  2Ar 

.  ^i,j+l  ~l~  V’ij— 1 

Az2 

Velocity  Field 


^  ~  Ti  2Az 

^  r-i  2Ar 


(4.20) 

(4.21) 

(4.22) 


Assuming  all  quantities  known  at  a  given  time  level,  the  solution  is  ad¬ 
vanced  by  one  time  step  as  summarized  below: 

1.  Equation  (4.9)  is  first  integrated  to  yield  new  values  of  the  vorticity 

everywhere  except  at  solid  boundaries. 

2.  Equation  (4.20)  is  inverted  using  Gauss-Seidel  iterations,  thus  yielding 

the  streamfunction  distribution  at  the  new  level. 

3.  The  no-slip  boundary  condition  is  imposed  using  Thom’s  approximation, 

giving  solid  wall  vorticity  boundary  conditions. 

4.  The  velocity  distribution  is  determined  using  Eqs.  (4.21  -  4.22),  and  the 

viscous  dissipation  function  at  the  new  time  level  is  computed. 

5.  The  temperature  distribution  at  the  new  time  level  is  obtained  by  inte¬ 

grating  the  energy  equation  (4.10)  and  imposing  boundary  conditions. 

4.2.1  Results  and  Discussion 

The  direct  simulation  code  is  first  applied  to  predict  low-speed  water  injec¬ 
tion  into  the  nozzle.  The  latter  is  discretized  on  a  rectangular  grid  with 
Nf  =  201  points  in  the  radial  direction  and  Ng  =  801  in  the  streamwise  di¬ 
rection.  A  20m/ s  injection  velocity  is  assumed,  corresponding  to  a  Reynolds 
number  Re  =  41, 626.  The  Prandtl  number  Pr  =  6.616,  and  the  Eckert  num¬ 
ber  Ec  =  0.000324  based  on  an  inlet  temperature  To  =  295°A.  The  inlet 
and  initial  velocity  and  temperature  profiles  are  adapted  from  ID  analytical 


expressions  of  Section  3,  assuming  an  inlet/initial  boundary  layer  thickness 
<^m  =  50.8fj,m.  A  10ms  injection  duration  is  simulated  using  an  integration 
time  step  At  =  0.002/is. 

Results  of  the  computations  are  shown  in  Figs.  4. 1-4.3,  which  respectively 
depict  streamwise  velocity  profiles  at  different  injection  times,  streamwise 
velocity  profiles  at  different  downstream  cross-sections,  and  temperature  dis¬ 
tributions  at  different  streamwise  planes.  At  early  injection  times,  the  figures 
illustrate  the  diffusive  growth  of  the  viscous  and  thermal  boundary  layers. 
For  small  downstream  locations,  z  <  10mm,  the  thermal  and  viscous  bound¬ 
ary  layers  tend  rapidly  towards  a  steady-state  value,  within  less  than  2ms 
after  the  start  of  the  injection.  Even  at  large  downstream  locations,  steady 
state  value  are  reached  by  the  end  of  the  injection  period.  Examination  of 
the  temperature  profiles  indicates  that,  for  the  present  injection  parameters, 
shear  heating  mechanisms  are  extremely  small.  Peak  temperature  values, 
which  are  recorded  near  the  nozzle  exit,  are  only  0.18"^'  higher  than  the  in¬ 
let  temperature.  It  is  interesting  to  note,  however,  that  most  of  the  heating 
actually  occurs  during  the  early  stage  of  injection,  during  which  the  bound¬ 
ary  layer  is  thin.  Furthermore,  comparison  of  the  temperature  profiles  at 
different  streamwise  locations  shows  that  the  wall  temperature  rises  rapidly 
at  small  downstream  locations,  where  the  boundary  layer  remains  thin,  and 
at  significantly  smaller  rates  as  we  move  downstream.  As  discussed  in  further 
detail  below  these  trends  will  also  be  observed  even  as  the  injection  speed  is 
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considerably  increased  and  shear  heating  effects  become  significantly  more 
pronounced. 

Unfortunately,  one  major  disadvantage  in  the  application  of  the  direct 
simulation  approach  is  the  high  computational  overhead  associated  with  the 
simulation.  For  the  low  injection  velocity  considered  above,  a  10ms  injec¬ 
tion  period  required  in  excess  of  two  weeks  to  complete  on  an  IBM  R6000- 
350  workstation.  For  higher  injection  velocities,  which  would  require  sig¬ 
nificantly  finer  resolutions,  the  necessary  computational  overhead  would  be 
prohibitively  high.  This  would  be  the  case  even  if  more  efficient  versions  of 
the  computations  were  implemented,  such  as  fast  solvers,  implicit  integration 
schemes,  or  more  accurate  spatial  discretizations. 

In  order  to  overcome  this  computational  difficulty,  an  alternative  ap¬ 
proach  is  adopted  which  is  based  on  the  parabolized  approximation  of  the 
equations  of  motion.  This  results  in  order-of-magnitude  savings  in  the  in¬ 
version  of  the  elliptic  streamfunction  operator,  which  dominates  the  over¬ 
all  computational  cost.  Further  reduction  of  CPU  requirements  is  sought 
through  the  implementation  of  stretched  computational  grids,  which  concen¬ 
trate  mesh  points  in  the  neighborhood  of  thin  computational  boundaries. 
The  construction  and  implementation  of  several  numerical  schemes  for  the 
simulation  of  the  parabolized  equations  of  motion  are  discussed  in  Section 
4.3  below. 
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t=1.0  ms 


t=4.0  ms 


Figure  4.1:  Laminar  boundary  development  in  an  axisymmetric  nozzle.  Wa¬ 
ter  is  injected  at  20m/s;  the  nozzle  is  20mm  long  and  has  2mm  inner  radius. 
The  plots  show  normalized  streamwise  velocity  profiles  at  different  times  fol¬ 
lowing  injection.  In  each  gragh,  profiles  at  different  downstream  locations 
are  plotted. 
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Figure  4.2:  Laminar  boundary  layer  development  in  an  axisymmetric  nozzle. 
Water  is  injected  at  2^ml s;  the  nozzle  is  20mm  long  and  2mm  inner  radius. 
The  plots  show  normalized  streamwise  velocity  profiles  at  different  stream- 
wise  locations.  In  each  gragh,  profiles  at  different  times  following  injection 
are  plotted. 
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Figure  4.3:  Thermal  boundary  layer  development  in  an  axisymmetric  nozzle. 
Water  is  injected  at  20m/s;  the  nozzle  is  20mm  long  and  has  2mm  inner 
radius.  The  plots  show  normalized  temperature  profiles  at  different  stream- 
wise  locations.  In  each  graph,  profiles  at  different  times  following  injection 
are  plotted. 


4.3  Parabolized  Approximations 


The  parabolized  approximation  is  motivated  by  the  fact  that  in  thin  bound- 
ary  layer  flows,  such  as  those  considered  here,  radial  diffusion  fluxes  dominate 
their  streamwise  counterparts.  Thus,  the  approximation  is  based  on  drop¬ 
ping  streamwise  diffusion  terms  from  the  governing  equations,  resulting  in 
the  system  summarized  in  Appendix  A. 


Numerical  simulation  of  the  parabolized  equations  of  motion  is  performed 
using  a  stretched  grid  technique.  To  this  end,  the  computational  (r,  z)  plane 
is  stretched  in  the  radial  direction  only,  i.e.  using  a  transformation  of  the 
form:  (r,z)  — >  (^,z).  In  all  the  computations  performed  in  this  study,  the 
transformation: 


^  ^  exp(a(l  -  ^))  -  1 

exp(o)  —  1 


(4.23) 


is  used.  Note  that  this  transformation  has  a  free  parameter,  a,  which  may 
be  adjusted  depending  on  the  inlet  boundary  layer  thickness  and  desired 
numerical  resolution.  Once  a  value  of  a  is  selected,  a  standard  rectangular 
finite-difference  grid  is  used  in  the  discretization  of  the  (^,  z)  domain.  The 
equations  of  motion  in  the  transformed  (^,  z)  plane  are  given  in  Appendix  A, 
which  also  discusses  the  numerical  simulation  of  this  equation  system. 


4.4  Results  and  Discussion 

4.4.1  Validity  of  the  parabolized  approximation 

Applications  of  the  parabolized  simulation  codes  starts  with  an  examination 
of  the  validity  of  the  approximation.  To  this  end,  steady  and  unsteady  sim¬ 
ulation  codes  are  tested  against  each  other  and  against  the  predictions  of 
direct  simulations  schemes.  The  same  injection  parameters  selection  in  Sec¬ 
tion  4.3  are  used. 

Comparison  of  wall  temperature  and  streamwise  velocity  predictions  ob¬ 
tained  using  the  steady  parabolized  equations  and  the  unsteady  parabolized 
equations  at  large  injection  times  {t  >  10ms)  reaveal  nearly  identical  results. 
Thus,  we  omit  discussion  of  the  unsteady  parabolized  equations,  and  focus 
on  contrasting  predictions  of  the  steady  parabolized  code  and  the  direct  sim¬ 
ulation  scheme. 

Fig.  4.4  shows  the  streamwise  growth  of  the  momentum  boundary  layer 
as  predicted  by  the  steady  parabolized  simulation  scheme.  This  simulations 
are  performed  on  an  unstretched  square  mesh  having  Nr  =  401  grid  points 
in  the  radial  direction  and  Nz  =  8001  grid  points  in  the  streamwise  direc¬ 
tion.  Comparison  of  these  results  with  those  of  Figs.  4.1-4.2  reveals  a  very 
favorable  agreement  between  both  the  approaches  at  large  times  following 
the  start  of  injection.  This  agreement  is  not  surprising  since,  as  mentioned 


in  Section  4.3,  the  boundary  layer  thickness  has  essentially  reached  its  steady 
value  at  the  late  stages  of  the  simulation. 

Further  comparison  between  the  two  prediction  schemes  is  presented  in 
Fig.  4.5,  which  depicts  streamwise  velocity  profiles  at  z  =  19.925mm  down¬ 
stream  of  the  nozzle  inlet.  The  profiles  are  drawn  using  results  of  the  steady 
parabolized  approximation  and  the  unsteady  direct  simulation  at  t  =  10ms 
following  injection.  Again,  an  excellent  agreement  between  the  two  ap¬ 
proaches  is  observed,  showing  only  small  deviations  between  the  results.  It 
is  also  interesting  to  note  that  the  steady  parabolized  approximation  pre¬ 
dicts  a  slightly  smaller  boundary  layer  thickness  near  the  nozzle  exit.  Thus, 
the  parabolized  approximation  is  not  expected  to  overpredict  boundary  layer 
growth  nor  underpredict  shear-heating  effects. 


4.4.2  Further  Analysis  of  Modeling  approximation 

The  efficiency  of  the  parabolized  numerical  simulation  schemes  enables  de¬ 
tailed  study  of  the  free  parameters  which  are  part  of  our  modeling  approach. 
As  indicated  earlier  in  this  section,  the  initial  and  inlet  conditions  used  in  the 
simulations  assume  a  quasi  ID  flow,  so  that  an  inlet  momentum  boundary 
layer  thickness,  and  an  initial  temperature  distribution  must  be  provided, 
thus,  it  is  desirable  to  first  examine  the  impact  of  free  parameters  on  the 
results  of  the  simulations  before  applying  the  codes  in  a  predictive  manner. 

This  exercise  is  conducted  for  high-speed  injection  of  kerosine.  Specif¬ 
ically,  an  injection  velocity  U  =  365m/s  and  an  inlet  temperature  = 
293°K  are  assumed  in  all  the  calculations  discussed  in  this  section.  The 
physical  properties  of  the  liquid  are  assumed  to  be  constant,  and  approx¬ 
imated  by  their  value  at  the  inlet  temperature;  we  use  p  -  393kg 
Cp  =  2093J/kg°K,  k  =  O.U88W/m°K,  and  u  =  2.263  x  IQ-^m^/s.  These 
injection  characteristics  are  thus  characterized  by  the  following  values  of 
Reynolds,  Prandtl  and  Eckert  numbers:  Re  =  322,576,  Pr  =  25.6,  and 
Ec  =  0.217. 


The  effect  of  the  inlet  momentum  boundary  layer  thickness,  S^,  is  first  ex¬ 
amined.  Both  steady  and  unsteady  computations  are  conducted  for  different 
values  of  Results  are  summarized  in  Fig.  4.7,  which  shows  steady-state 
wall  temperature  distributions  for  all  cases.  A  flat  inlet  temperature  profile 
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Figure  4.4:  Streamwise  velocity  profiles  for  the  same  injection  parameters  of 
Fig.  4.1,  computed  using  steady  parabolized  simulations.  The  streamwise 
location  is  indicated. 
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Figure  4.5:  Streamwise  velocity  profiles  at  2:  =  19.925mm  for  the  same  injec¬ 
tion  parameters  of  Fig.  4.2,  computed  using  steady  paravolized  simulations 
and  the  unsteady  direct  simulation  at  t  =  10ms. 


is  assumed,  and  the  initial  temperature  distribution  is  taken  to  be  uniform. 
Fig.  4.7  shows  that,  while  the  temperature  behavior  near  the  nozzle  inlet 
may  be  strongly  affected  by  the  value  of  Sin,  heating  characteristics  further 
downstream  only  show  a  weak  dependence  on  inlet  conditions.  In  partic¬ 
ular,  as  the  inlet  momentum  boundary  layer  thickness  decreases,  the  peak 
wall  temperature  reached  at  the  nozzle  exit  exhibits  small  variation  with 
Sin-  Specifically,  when  the  inlet  momentum  boundary  layer  thickness  be¬ 
comes  of  the  order  of  Ifim  or  smaller,  deviations  in  the  peak  temperature 
within  the  field  are  less  than  5°K.  This  deviation  is  a  very  small  fraction 
of  the  total  variation  of  the  wall  temperature,  which  increases  by  more  than 
U0°K  across  the  length  of  the  nozzle.  Consequently,  the  effect  of  the  inlet 
momentum  boundary  layer  thickness  may  be  safely  absorbed  by  consistently 
decreasing  5^  until  peak  temperature  predictions  are  essentially  independent 
of  selected  values. 

In  order  to  further  support  this  assessment,  the  predictions  of  steady  and 
unsteady  parabolized  approximations  are  tested  against  each  other.  This 
exercise  is  summarized  in  Fig.  4.8  which  shows  instantaneous  wall  tempera¬ 
ture  distributions  for  the  same  conditions  summarized  above,  starting  with 
an  inlet  momentum  boundary  layer  thickness  Sin  =  1. 29 fxm.  The  computed 
results  indicate  that,  at  early  stage,  a  near-uniform  heating  of  the  fluid  occurs 
except  near  the  nozzle  inlet  where  spatial  variations  of  the  boundary  layer 
thickness  are  important.  Thus,  during  early  stages  of  injection,  the  heating 


at  large  downstream  locations  occurs  in  a  quasi  one-dimensional  fashion,  as 
assumed  in  the  analysis  of  section  3. 

For  larger  times,  t  >  1.5  ms,  the  thermal  and  viscous  boundary  layers 
reach  their  steady  state  values  in  the  entire  nozzle.  Comparison  of  the  wall 
temperature  distribution  for  t  >  1.5  ms  with  corresponding  steady  state  pre¬ 
dictions  (Fig.  4.7)  reveal  nearly  identical  results.  An  isolated  comparison 
is  thus  omitted,  since  the  time  interval  required  to  reach  a  steady-state  is 
smaller  than  the  injection  duration  of  interest,  the  application  of  a  steady 
analysis  to  the  prediction  of  peak  temperatures  proves  sufficient.  However, 
the  application  of  unsteady  simulation  codes  is  still  performed  in  most  ap¬ 
plications  described  in  this  work,  primarily  as  an  additional  means  of  check¬ 
ing  the  computed  predictions.  This  constitutes  a  valid  approach  since,  as 
explained  in  Appendix  A,  steady  and  unsteady  simulation  codes  rely  on  dif¬ 
ferent  discretization  and  integration  methodologies. 

Finally,  the  impact  of  inlet  thermal  profile  and  initial  temperature  distri¬ 
bution  is  analyzed.  To  this  end,  predictions  based  on  a  flat  inlet  temperature 
profile  and  uniform  initial  temperature  distribution  are  contrasted  to  those 
obtained  using  the  quasi-lD  results  to  prescribe  the  inlet  profile  and  initial 
distribution.  The  comparison,  performed  using  steady  parabolized  approx¬ 
imation  with  5in  =  1.29/im,  is  shown  in  Fig  4.9.  The  figure  indicates  that, 
despite  significant  differences  at  small  downstream  locations,  steady  state 


predictions  of  the  peak  wall  temperature  near  the  nozzle  exit  are  weakly 
sensitive  to  inlet  conditions.  When  coupled  with  the  results  of  the  above 
analysis,  this  enables  us  to  conclude  that  peak  temperature  predictions  are 
essentially  unsensitive  to  both  thermal  and  viscous  inlet  conditions. 


4.4.3  Shear  Heating  of  LP  1846  during  High-Speed  In¬ 
jection 

In  this  section,  the  parabolized  simulation  schemes  are  applied  to  charac¬ 
terize  shear  induced  heating  during  high-speed  injection  of  liquid  monopro¬ 
pellant  LP  1846.  In  the  computations  of  this  section,  an  inlet  temperature 
To  =  298‘’ii"  is  considered  in  all  cases.  Furthermore,  the  physical  properties  of 
LP  1846  are  assumed  constant  and  approximated  by  their  value  at  the  injec¬ 
tion  temperature  (see  Section  5,  below).  Specifically,  we  use  p  =  1, 400A:p/m^, 
Cp  =  2, 300J/kg°K,  k  =  0.15W/m°K,  and  u  =  4.988  x 

The  impact  of  injection  velocity,  is  studied  by  considering  four  different 
values:  U  =  lOOm/s,  200m/s,  300m/s,  400m/s.  The  corresponding  injec¬ 
tion  characteristics  are  respectively  characterized  by  the  following  Reynolds- 
Eckert  number  pairs:  {Re  =  40096,  Ec  =  0.014},  [Re  =  80191,  Ec  =  0.058}, 
[Re  =  120287,  Ec  =  0.13},  and  {Re  =  160384,  Ec  =  0.23}.  Since  constant 
property  models  are  applied,  all  cases  are  characterized  by  the  same  Prandtl 
number,  Pr  =  109.4. 

Results  of  steady  computations  are  summarized  in  Fig.  4.10,  which  shows 
the  wall  temperature  distribution  for  all  four  cases.  As  in  Section  4.4.2,  the 
results  are  compared  to  predictions  of  the  unsteady  computations  in  order 
to  check  their  validity.  An  illustration  of  this  exercise  is  given  in  Fig.  4.11, 
which  shows  instantaneous  wall  temperature  distributions  for  an  injection 


50 


Figure  4.6:  Wall  temperature  distributions  during  high-speed  injection  of 
kerosine  with  U  =  365m/ s  and  To  —  29S°K.  The  computations  are  per¬ 
formed  using  the  steady  parabolized  approximation  on  a  stretched  grid  woth 
Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the  stream- 
wise  direction.  The  grid  stretching  parameter  a  =  6.5,  amd  a  flat  inlet  tem¬ 
perature  profile  is  imposed.  The  inlet  momentum  boundary  layer  thickness 
Sin  is  indicated. 
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Figure  4.7:  Instantaneous  wall  temperature  distributions  during  high-speed 
injection  of  kerosine  with  U  =  365m/s  and  To  =  293°K.  The  computations 
are  performed  using  the  unsteady  parabolized  approximation  on  a  stretched 
grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the  radial 
direction  and  =  8001  points  in  the  streamwise  direction.  A  flat  inlet 
temperature  profile  is  imposed,  while  the  inlet  viscous  layer  has  =  1.29//m. 
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Figure  4.8:  Steady  state  wall  temperature  distributions  for  high-speed  injec¬ 
tion  of  kerosine  with  U  =  365m/ s  and  Tp  —  293°K.  The  computations  are 
performed  using  the  steady  parabolized  approximation  on  a  stretched  grid 
with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction 
and  =  8001  points  in  the  streamwise  direction.  The  inlet  viscous  layer 
has  Sin  =  1.29fj,m.  Results  are  shown  for  both  flat  (dashed  line)  and  non- 
uniform  (solid  line)  profiles.  In  the  latter  case,  the  inlet  temperature  profile 
is  adapted  from  the  analysis  of  section  3. 


velocity  of  300m/s.  For  large  times  following  injection,  t  >  1.33ms,  the  wall 
temperature  distribution  in  the  unsteady  computation  coincides  with  steady 
state  prediction.  Thus,  a  steady  field  is  rapidly  reached,  well-before  the  end 
of  injection  period. 

The  computations  shown  in  Fig.  4.10  also  reflect  the  trends  of  quasi  ID 
analysis.  In  particular,  the  predictions  also  exhibit  a  quadratic  dependence 
on  the  injection  velocity.  Significant  heating  of  the  mixture,  with  tempera¬ 
ture  increases  greater  than  a  im°K,  are  predicted  when  the  injection  velocity 
exceeds  200m/ s.  Thus,  for  these  injection  scenarios,  premature  ignition  of 
the  mixture  due  to  severe  shear  heating  is  likely  to  occur. 

It  is  interesting  to  note  that  shear  heating  characteristics  for  LP  1846 
are  of  the  same  order  as  that  computed  for  kerosene.  This  similarity  may 
be  observed  by  comparing  Figs.  4.10  and  4.7.  The  comparison  reveals  that 
the  peak  wall  temperature  for  a  300m/s  injection  of  LP  1846  is  close  to  that 
achieved  for  kerosene  injection  at  365m/ s.  Thus,  shear  heating  mechanisms 
of  LP1846  are  more  pronounced  than  those  of  kerosene. 

It  is  also  interesting  to  note  that  the  300m/s  injection  of  LP  1846  is  char¬ 
acterized  by  a  smaller  Eckert  number  than  kerosene  injection  at  365m/s,  and 
that  the  inlet  temperatures  in  both  cases  are  nearly  identical.  On  the  other 
hand,  the  Prandtl  number  of  kerosene  is  significantly  smaller  than  that  of  LP 


1846,  which  may  explain  the  more  pronounced  heating  effects  observed  in  the 
latter  case.  Further  discussion  of  the  role  of  Prandtl,  Eckert  and  Reynolds 
numbers  is  provided  in  Section  5  below. 


4.4.4  Effect  of  Wall  Heat  Transfer 

Finally,  the  impact  of  wall  heat  transfer  on  peak  temperature  prediction  is 
investigated.  We  use  a  simplified  model  in  which  wall  heat  transfer  is  taken 
into  account  through  a  convection  heat  transfer  coefiicient.  Thus,  the  sim¬ 
plified  model  ignores  the  thermal  resistance  of  finite  thickness  nozzle  walls, 
and  also  ignores  the  associated  heat  storage  capacity  which  may  play  an 
important  role  during  the  flow  transient.  Consequently,  the  model  assumes 
that  the  nozzle  walls  are  extremely  thin  and  that  the  limiting  heat  transfer 
mechanism  from  the  liquid  to  its  surroundings  is  due  to  convection  from  the 
nozzle’s  outer  boundaries. 

In  the  computations,  convection  heat  transfer  from  the  nozzle  walls  is 
incorporated  by  modifying  thermal  boundary  condition  at  the  nozzle  radius. 
The  modified  boundary  conditions,  which  expresses  the  continuity  of  the  heat 
flux,  is  written  as: 
dT 

~^~^\r=Ro  —  ^(^|r=iio  ~  Too)  (4.24) 

where  k  is  the  thermal  conductivity  of  the  mixture,  h  is  the  heat  transfer 
coefficient  and  is  the  far-field  temperature  of  the  surrounding  fluid.  In  all 


computations,  we  assume  that  the  far-field  temperature  of  the  surrounding 
fluid  coincides  with  the  inlet  mixture  temperature,  i.e.  Too  =  Tg.  Substituting 
this  assumption  into  Eq.  (4.24)  and  normalizing  the  resulting  expression,  we 
get: 

dT 

=  1)  =  Nu(T{r  =  1)  -  1)  (4.25) 

where 


Nu  =  ^ 
k 


(4.26) 


is  the  Nusselt  number  based  on  the  mixture’s  thermal  conductivity.  Note  that 
insulated  wall  conditions  may  also  be  simulated  simply  by  setting  Nu  =  0, 
in  which  case  Eq.  (4.25)  reduces  to  the  homogeneous  Neumann  Boundary 
condition. 


High  speed  LP  1846  injection  experiments  are  conducted  for  diflferent 
values  of  the  heat  transfer  coefficient.  We  select  three  characteristic  values 
of  the  heat  transfer  coefficient,  h  =  20W/rn^.^K,  h  =  lOOW/ni^.^K,  and 
h  =  b00W/m^.°K.  These  values  are  respectively  representative  of  free  con¬ 
vection  conditions  in  air,  forced  air  cooling  at  low  speed,  and  forced  liquid 
cooling  at  moderate  speed.  The  corresponding  Nusselt  numbers  defined  by 
Eq.  (4.26)  are  Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67  respectively. 

Computed  results  are  shown  in  Figs.  4.11-4.13,  which  respectively  show 
peak  temperature  distributions  for  LP  1846  injection  at  U  =  100m/ s,  U  = 


200m/s„  and  U  =  300m/s.  In  each  plot,  curves  are  drawn  for  insulated  wall 
conditions  {Nu  =  0)  and  for  three  Nusselt  number  values  specified  above. 
Examination  of  these  predictions  reveals  that: 

1.  For  very  high  injection  speed  {U  =  300m/s)  wall  heat  transfer  does  not 

significantly  reduce  peak  temperature  predictions  even  at  high  Nus¬ 
selt  number.  Thus,  conventional  nozzle  wall  cooling  means  may  not 
constitute  an  effective  means  of  minimizing  the  likelihood  of  mixture 
preignition. 

2.  For  low  Nusselt  number,  Nu  =  267,  the  computed  peak  temperatures 

are  very  close  to  those  obtained  assuming  adiabatic  wall  conditions. 
Accordingly,  natural  heat  convection  from  small  diameter  nozzles  is 
not  expected  to  significantly  affect  peak  temperature  predictions. 

3.  Large  values  of  the  heat  transfer  coefficient  may  appreciably  reduce  peak 

temperatures  in  thin  walled  nozzles  whenever  the  injection  velocity  is 
not  extremely  high.  In  these  situations,  forced  cooling  techniques  may 
be  especially  tailored  in  order  to  effectively  minimize  the  risk  of  mixture 
ignition. 

It  is  finally  emphasized  that,  when  wall  heat  transfer  is  accounted  for,  the 
peak  temperature  achieved  at  a  given  streamwise  location  may  not  always 
coincide  with  the  wall  temperature.  Generally,  as  heat  losses  increase,  the 
maximum  temperature  location  moves  away  from  the  wall  into  the  thermal 
boundary  layer.  This  effect  is  illustrated  in  Figs  4.14-4.15,  which  respectively 


show  radial  temperature  profiles  at  the  nozzle  exit  for  injection  velocities 
U  =  100m/ s  and  300m/ s.  Consequently,  unlike  insulated  wall  conditions, 
the  peak  temperature  distributions  plotted  in  Figs.  4.11-4.13  do  not  always 
correspond  to  wall  temperature  distributions. 


Figure  4.9:  Steady  state  wall  temperature  distributions  for  high-speed  injec¬ 
tion  of  LP  1846  with  To  =  298°/i£r  at  four  different  velocities,  U  =  lOOm/s, 
U  —  200m/s,  U  —  300m/s,  and  U  =  400m/s.  The  computations  are  per¬ 
formed  using  the  steady  parabolized  approximation  on  a  stretched  grid  with 
a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and 
Nz  =  8001  points  in  the  streamwise  direction. 
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Figure  4.10:  Instantaneous  wall  temperature  distribution  for  high-speed  in¬ 
jection  of  LP  1846  with  =  298° K  and  U  =  SOOm/s.  The  computations  are 
performed  using  the  unsteady  parabolized  approximation  on  a  stretched  grid 
with  a  stretching  parameter  a  =  6,  Nr  =  201  points  in  the  radial  direction 
and  Nz  =  801  points  in  the  streamwise  direction.  The  integration  time  step 
At  =  4  X  10-1 
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Figure  4.11;  Steady  state  peak  temperature  distributions  for  high-speed  in¬ 
jection  of  LP  1648  with  To  =  298"Jf,  U  =  lOOm/s,  and  four  different  Nusselt 
number  Nu  =  0,  ATu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  com¬ 
putations  are  performed  using  the  steady  parabolized  approximation  on  a 
stretched  grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the 
radial  direction  and  —  8001  points  in  the  streamwise  direction. 
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Figure  4.12;  Steady  state  peak  temperature  distributions  for  high-speed  in¬ 
jection  of  LP  1648  with  To  =  298"K,  U  =  200m/s,  and  four  different  Nusselt 
number  Nu  =  0,  Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  com¬ 
putations  are  performed  using  the  steady  parabolized  approximation  on  a 
stretched  grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the 
radial  direction  and  =  8001  points  in  the  streamwise  direction. 
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Figure  4.13:  Steady  state  peak  temperature  distributions  for  high-speed  in¬ 
jection  of  LP  1648  with  To  =  298° A',  U  =  300m/s,  and  four  different  Nusselt 
number  Nu  =  0,  Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  com¬ 
putations  are  performed  using  the  steady  parabolized  approximation  on  a 
stretched  grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the 
radial  direction  and  =  8001  points  in  the  streamwise  direction. 
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z  =  2.0  cm 


Figure  4.14:  Steady  state  radial  temperature  profile  at  the  nozzle  exit  for 
high-speed  injection  of  LP  1648  with  To  =  298°/^,  U  =  lOOm/s,  and  four 
different  Nusselt  number  Nu  =  0,  iVu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67. 
The  computations  are  performed  using  the  steady  parabolized  approximation 
on  a  stretched  grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in 
the  radial  direction  and  =  8001  points  in  the  streamwise  direction. 
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z  =  2.0  cm 


Figure  4.15:  Steady  state  radial  temperature  profile  at  the  nozzle  exit  for 
high-speed  injection  of  LP  1648  with  To  =  298° K,  U  =  300m/ s,  and  four 
different  Nusselt  number  Nu  =  0,Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67. 
The  computations  are  performed  using  the  steady  parabolized  approximation 
on  a  stretched  grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in 
the  radial  direction  and  =  8001  points  in  the  streamwise  direction. 


Chapter  5 

Variable  Mixture  Properties 


In  this  section,  the  constant-properties  assumption  is  relaxed  and  shear  heat¬ 
ing  mechanisms  are  examined  for  a  mixture  with  temperature-dependent  vis¬ 
cosity  and  thermal  conductivity.  The  nozzle  geometry  and  injection  regime 
of  interest  are  identical  to  those  considered  in  the  previous  section.  However, 
the  accommodation  of  a  variable  property  mixtures  necessitates  a  signifi¬ 
cant  change  both  in  the  governing  equations  (given  in  section  5.1)  and  in 
the  corresponding  numerical  simulation  schemes.  The  latter  are  discussed 
in  Appendix  B,  which  summarizes  both  steady  and  unsteady  schemes  for 
the  simulation  of  the  parabolized  equations  where  the  temperature  depen¬ 
dence  of  the  viscosity  and  thermal  conductivity  may  not  be  exactly  known. 
In  section  5.3,  the  computational  codes  are  applied  to  predict  shear-induced 
heating  of  liquid  monopropellant  LP  1846  for  different  injection  speeds,  inlet 
temperatures,  and  wall  cooling  conditions. 
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5 . 1  Formulation 


As  in  Section  4,  we  nssume  thnt  the  mixture  is  an  incompressible  liquid, 
and  focus  on  the  same  injection  parameters  and  nozzle  geometries  consid¬ 
ered  there.  However,  the  viscosity  and  thermal  conductivity  of  the  mix¬ 
ture  are  now  allowed  to  vary  with  temperature.  This  generalization  does 
not  necessitate  any  changes  in  the  basic  modeling  of  the  flow  domain  or 
boundary  conditions.  Furthermore,  all  kinematical  relationships  are  left  un¬ 
affected.  Specifically,  vorticity-streamfunction  relationships  and  the  viscous 
dissipation  function  definition  introduced  in  previous  sections  continue  to 
hold.  Thus,  we  limit  our  discussion  to  a  simple  statement  of  the  vorticity 
transport  and  energy  conservation  equations,  whose  form  changes  in  order 
to  accommodate  the  desired  extension. 


Using  the  same  normalization  conventions  introduced  in  Section  2,  and 
introducing  the  following  normalized  viscosity  and  thermal  conductivity  def¬ 
initions. 


k* 


=  £21 

Hfo) 

-  Hf) 

Hfo) 


(5.1) 

(5.2) 


where  ~  denotes  dimensional  quantities,  the  vorticity  transport  and  en¬ 
ergy  conservation  equations  are  respectively  expressed  as: 
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Here,  the  Reynolds  and  Prandtl  numbers  are  based  on  the  viscosity  and 


thermal  conductivity  measured  at  the  inlet  mixture  temperature,  i.e. 


Re  = 
Pr  = 


tVR 
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5.2  Modeling  Strategies 


The  applications  discussed  in  this  section  are  motivated  by  a  desire  to  char¬ 
acterize  the  likelihood  of  preignition  due  to  shear-induced  heating  of  liquid 
monopropellants  during  high-speed  injection;  LP  1846  will  be  specifically  se¬ 
lected  in  all  applications.  For  this  mixture,  the  dependence  of  the  dynamic 
viscosity  on  temperature  may  be  expressed  as: 


fl  =  C  exp 


D 


.T-Tr 


ref  ^ 


(5.7) 


where  pi  is  the  viscosity  measured  in  centipoise(cp),  f  is  the  temperature  mea¬ 
sured  in  Kelvin,  Tref  =  164°K  is  a  reference  temperature  while  C  =  0.16773 


cp  and  D  —  b02.52°K  are  dimensionalconstants.  Unfortunately,  the  depen¬ 
dence  of  the  thermal  conductivity  of  LP-1846  on  the  prevailing  temperature 
is  not  known,  and  only  a  single  value  k{f  =  298° K)  =  0.15W/m°K  is  re¬ 
ported  in  the  literature. 

In  order  to  tackle  the  uncertainty  regarding  temperature-related  vari¬ 
ations  of  the  thermal  conductivity,  two  modeling  approaches  are  adopted 
in  the  following  computations.  The  first  modeling  approach  is  based  on 
the  observation  that  in  most  liquids,  variations  of  the  thermal  conductiv¬ 
ity  are  much  smaller  than  those  of  the  dynamic  viscosity.  Specifically,  for 
most  liquids,  the  dynamic  viscosity  decreases  rapidly  with  increasing  tem¬ 
perature,  with  corresponding  small  variation  in  the  thermal  conductivity. 
Accordingly,  the  Prandtl  number  is  expected  to  drop  appreciably  with  in¬ 
creasing  temperature.  Thus,  the  first  model  calls  for  treating  the  thermal 
conductivity  as  temperature-independent,  and  using  the  reported  value  for 
a  mixture  at  298" A”.  Since  both  the  density  and  heat  capacity  of  LP-1846 
(cp  =  280QJ/kg°K)  vary  slightly  with  temperature,  the  thermal  diffusivity  is 
consequently  constant  and  the  first  model  is  called  the  uniform  Peclet  num¬ 
ber  model. 

The  second  modeling  approach  is  motivated  by  the  theoretical  findings 
of  Section  3  and  the  simulation  of  section  4,  which  indicate  that,  other  pa¬ 
rameters  being  equal,  shear  heating  effects  tend  to  be  significantly  more  pro- 


nounced  for  mixtures  having  higher  Prandtl  number.  Thus,  the  second  model 
conservatively  assumes  that  the  thermal  conductivity  (and  consequently  the 
thermal  diffusivity)  admits  a  temperature  dependence  which  is  similar  to 
that  of  the  dynamic  viscosity.  Thus,  the  second  model  is  called  the  uniform 
Prandtl  number  model.  It  is  expected  to  yield  more  conservative  estimated 
peak  temperatures  since  it  ignores  the  potential  decrease  in  the  Prantle  num¬ 
ber  with  increasing  temperature.  Both  models  are  evaluated  in  the  simula¬ 
tions  discussed  below. 


5.3  Application  of  LP-1846 

Simulation  schemes  accommodating  variable-property  mixtures  are  applied 
to  re-examine  the  predictions  shown  in  Section  4.4.3.  High-speed  injection  of 
LP-1846  is  once  again  considered,  with  injection  velocities,  U  =  lOOm/s, 
U  —  200m/s,  U  =  30077i/s,  and  U  =  400m/s.  In  all  cases,  an  injec¬ 
tion  temperature  To  =  29S°K  is  assumed.  These  injection  experiments  are 
respectively  characterized  by  the  following  Reynolds-Eckert  number  pairs: 
{Re  =  40096,  Ec  =  0.014),  {Re  =  80109,  Ec  =  0.058),  {Re  =  120287, 
Ec  =  0.13),  and  {Re  =  160384,  Ec  =  0.23),  and  by  a  Prandtl  Pr  =  109.4. 
All  dimensionless  groups  are  based  on  properties  evaluated  at  the  injection 
temperature.  Unless  otherwise  stated,  adiabatic  wall  conditions  are  assumed. 


Results  of  the  computations  are  summarized  in  Figs.  5.1  and  5.2,  which 


respectively  show  steady  state  wall  temperature  distributions  obtained  using 
the  uniform  Peclet  and  Prandtl  number  models.  As  before,  the  validity  of 
these  predictions  are  checked  against  those  of  unsteady  computations.  Since 
a  similar  agreement  to  the  observed  in  Section  4.4  is  again  observed,  results 
of  this  exercise  are  omitted. 

Comparison  of  the  present  results  with  those  obtained  in  Section  4.4.3 
are  summarized  as  follows: 

(1)  For  moderate  injection  speeds,  U  <  200m/s,  the  predictions  of  the  uni¬ 

form  Peclet  number  model  are  very  close  for  those  obtained  assuming 
constant  properties.  At  higher  injection  speeds,  the  larger  temperature 
variations  induced  by  intense  shear  heating  of  the  mixture  cause  a  sig¬ 
nificant  deviation  between  the  predictions.  The  uniform  Peclet  number 
model  predicts  lower  steady  state  peak  temperatures  than  the  corre¬ 
sponding  constant  property  simulation.  The  nature  of  this  deviation 
is  not  surprising  since  the  viscosity  decreases  with  increasing  tempera¬ 
ture,  resulting  in  a  drop  in  the  Prandtl  number.  Thus,  by  neglecting  the 
temperature  dependence  of  the  viscosity,  the  constant  property  model 
yields  more  conservative  estimates  of  shear  induced  heating  effects. 

(2)  When  the  Prandtl  number  is  artificially  held  constant,  an  unrealistically 

large  shear  heating  of  the  mixture  is  predicted.  Figure  5.2  shows  that 
the  uniform  Prandtl  number  model  yields  peak  temperature  which, 
for  large  injection  speeds,  may  be  twice  as  large  as  those  predicted 


by  the  other  two  models.  These  temperature  increases  are  deemed 
unrealistic  since,  as  previously  mentioned,  the  thermal  conductivity 
drops  slightly  with  increasing  temperature  while  the  viscosity  decreases 
more  substantially.  Thus,  by  prescribing  a  temperature  dependence  of 
the  thermal  conductivity  which  is  similar  to  that  of  the  viscosity,  the 
uniform  Prandtl  model  generally  yields  overly  conservative  estimates 
of  shear  heating  effects  and  should  not,  therefore,  be  relied  upon  as  a 
reliable  predictive  tool. 

Following  the  above  discussion,  only  the  uniform  Peclet  number  model  is 
used  in  the  computations  presented  below. 

5.3.1  Effect  of  Wall  Heat  Transfer 

The  impact  of  wall  heat  transfer  is  analyzed  in  a  similar  fashion  to  that 
adopted  in  constant  property  simulations.  We  use  the  uniform  Peclet  num¬ 
ber  model  and  apply  steady  parabolized  approximation  to  predict  shear  heat¬ 
ing  with  heat  transfer  conditions  characterized  by  the  following  coefficients: 
h  =  2m/vn?°K,  h  =  \mW/w?°K,  and  h  =  Note  that, 

since  a  uniform  Peclet  number  model  is  used,  the  thermal  conductivity  is 
taken  to  be  constant.  Consequently,  the  Nusselt  number  definition  given  in 
Section  4.3  also  holds,  and  modification  of  the  adiabatic  wall  condition  is 
performed  in  an  identical  manner  to  that  described  there. 


Results  of  the  computations  are  summarized  in  Figs.  5.3  —  5.5,  which 
show  that  peak  temperature  distributions  for  high  speed  LP  1846  injection 
with  To  =  298°iir,  and  C/  =  lOOm/s,  U  -  200m/ s,  and  U  =  300m/s.  Ex¬ 
amination  of  these  predictions  confirm  earlier  expectations  regarding  both 
the  role  of  heat  transfer  and  the  eff'ect  of  variable  properties.  Specifically, 
all  the  trends  established  using  the  constant-property  model  are  once  again 
observed.  Moreover,  comparison  of  the  results  of  the  uniform  Peclet  number 
and  constant-property  model  are  also  in  agreement  with  adiabatic  wall  pre¬ 
dictions.  In  particular,  when  the  injection  speed  is  low,  results  of  the  uniform 
Peclet  number  model  are  very  close  to  those  obtained  using  constant-property 
model.  For  higher  injection  speeds,  temperature  variations  are  significantly 
more  pronounced  and  the  uniform  Peclet  number  model  yields  smaller  peak 
temperature  predictions  than  those  obtained  using  constant-property  simu¬ 
lations. 


5.3.2  Effect  of  injection  temperature 

Finally,  the  effect  of  inlet  temperature  is  investigated.  We  consider  three 
injection  speeds,  U  =  lOOm/s,  U  =  200m/s,  and  U  =  300m/s,  and  assume 
adiabatic  nozzle  wall  conditions.  Results  of  steady  parabolized  approxima¬ 
tions  are  shown  in  Figs.  5.6  —  5.8,  which  respectively  show  wall  temperature 
distributions  for  three  inlet  temperatures.  To  =  278°K,  To  =  298'’^',  and 
To  =  318°K. 


For  low  injection  speed,  U  =  lOOm/s,  the  effects  of  shear-induced  heating 
are  essentially  similar  for  all  inlet  temperatures.  For  these  injection  charac¬ 
teristics,  moderate  heating  in  the  thermal  boundary  layer  occurs,  and  the 
wall  temperature  distributions  for  diflferent  cases  appear  to  be  shifted  ver¬ 
tically  as  the  inlet  temperature  is  varied.  This  result  is  not  surprising,  and 
is  in  agreement  with  previous  results,  since  the  temperature  distributions  - 
and  consequently  the  viscosity  -  do  not  exhibit  large  variations. 

At  higher  injection  speeds,  U  =  200m/s,  variable  viscosity  effects  start 
becoming  more  pronounced.  Figure  5.7  indicates  that  shear  heating  of  the 
mixture  is  more  substantial  as  the  inlet  temperature  is  decreased.  Note  that, 
for  lower  inlet  temperatures,  the  inlet  viscosity  of  the  mixture  is  higher. 
Therefore,  shear  stresses  and  viscous  dissipation  are  also  higher;  this  results 
in  larger  wall  temperature  increases.  This  trend  can  also  be  interpreted  in 
terms  of  the  expectation,  established  earlier  in  the  context  of  quasi  ID  flow 
and  constant-property  simulations,  that  shear  heating  effects  are  pronounced 
for  mixtures  with  higher  Prandtl  number.  The  present  results  are  consistent 
with  this  trend,  since  the  inlet  Prandtl  number  increases  with  decreasing 
temperature. 

Another  important  observation  is  that,  at  high  injection  velocities,  all 
similarity  between  wall  temperature  distributions  is  lost  as  the  inlet  temper- 


ature  developed  at  low  injection  temperatures  may  exceed  that  corresponding 
to  higher  inlet  temperature.  Further  examination  of  the  computations  indi¬ 
cates  that  these  trends  are  due  to  different  development  of  both  the  thermal 
and  viscous  boundary  layers.  Specifically,  as  illustrated  in  Figs.  5.9-5.10, 
when  intense  shear-heating  of  the  mixture  occurs,  both  the  structure  and 
spatial  evolution  of  the  boundary  layer  exhibit  significant  differences  as  the 
inlet  temperature  is  varied. 

The  observed  dependence  of  shear  heating  effects  on  the  injection  temper¬ 
ature  appears  to  pose  a  significant  challenge  to  the  evaluation  of  the  likelihood 
of  mixture  ignition.  However,  it  should  be  noted  that  this  dependence  admits 
a  consistent  trend,  namely  that  injection  characteristics  having  higher  inlet 
Prandtl  number  exhibit  higher  temperature  increases.  Thus,  conservative  es¬ 
timates  of  shear  heating  effects  can  still  be  obtained  by  focusing  the  analysis 
on  injection  conditions  which  are  characterized  by  the  highest  anticipated 
Prandtl  number. 


Figure  5.1:  Steady  state  wall  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  To  =  298° K  at  four  difference  velocities,  U  =  lOOm/s, 
U  =  200m/s,  U  =  SOOm/ s,  and  U  =  40077i/s.  The  uniform  Peclet  num¬ 
ber  model  is  adopted,  and  the  computations  are  performed  using  the  steady 
parabolized  approximation  on  a  stretched  grid  with  a  stretching  parameter 
a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the 
streamwise  direction. 
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Figure  5.2:  Steady  state  wall  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  Tq  =  298° K  at  four  difference  velocities,  U  =  lOOm/s, 
U  =  299m/ s,  U  =  300m/s,  and  U  =  400m/s.  The  uniform  Prandtl  num¬ 
ber  model  is  adopted,  and  the  computations  are  performed  using  the  steady 
parabolized  approximation  on  a  stretched  grid  with  a  stretching  parameter 
a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the 
streamwise  direction. 
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Figure  5.3:  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  %  =  298°K,  U  =  lOOm/s,  and  four  different  Nusselt 
number  Nu  =  0,  Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  uniform 
Peclet  number  model  is  adopted,  and  the  computations  are  performed  using 
the  steady  parabolized  approximation  on  a  stretched  grid  with  a  stretching 
parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001 
points  in  the  streamwise  direction. 
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Figure  5.4:  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  To  =  298"/f,  U  =  200m/s,  and  four  different  Nusselt 
number  Nu  =  0,  iVu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  uniform 
Peclet  number  model  is  adopted,  and  the  computations  are  performed  using 
the  steady  parabolized  approximation  on  a  stretched  grid  with  a  stretching 
parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001 
points  in  the  streamwise  direction. 
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Figure  5.5:  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  %  =  298°K,  U  =  300m/s,  and  four  different  Nusselt 
number  Nu  =  0,  Nu  =  0.267,  Nu  =  1.33,  and  Nu  =  6.67.  The  uniform 
Peclet  number  model  is  adopted,  and  the  computations  are  performed  using 
the  steady  parabolized  approximation  on  a  stretched  grid  with  a  stretching 
parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001 
points  in  the  streamwise  direction. 
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Figure  5.6;  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  U  =  100m/ s,  and  three  different  inlet  temperatures 
To  =  27S^K,  To  =  298‘’ii:,  and  To  =  S1S°K.  The  uniform  Peclet  num¬ 
ber  model  is  adopted,  and  the  computations  are  performed  using  the  steady 
parabolized  approximation  on  a  stretched  grid  with  a  stretching  parameter 
a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the 
streamwise  direction. 


z  (cm) 


Figure  5.7:  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  U  =  200m/s,  and  three  different  inlet  temperatures 
To  =  278" A',  To  =  298° K,  and  To  =  318° AT.  The  uniform  Peclet  num¬ 
ber  model  is  adopted,  and  the  computations  are  performed  using  the  steady 
parabolized  approximation  on  a  stretched  grid  with  a  stretching  parameter 
a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the 
streamwise  direction. 


Figure  5.8:  Steady  state  peak  temperature  distribution  for  high-speed  injec¬ 
tion  of  LP  1648  with  U  =  300m/s,  and  three  different  inlet  temperatures 
To  =  278° K,  To  =  298° K,  and  To  =  318° A!".  The  uniform  Peclet  num¬ 
ber  model  is  adopted,  and  the  computations  are  performed  using  the  steady 
parabolized  approximation  on  a  stretched  grid  with  a  stretching  parameter 
o  =  6.5,  Nr  =  401  points  in  the  radial  direction  and  =  8001  points  in  the 
streamwise  direction. 
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Figure  5.9;  Steady  state  streamwise  velocity  profiles  at  different  streamwise 
locations  for  high-speed  injection  of  LP  1846  with  U  =  300to/s,  and  two 
different  inlet  temperature  To  =  278°iC(dashed  line),  and  To  =  298° K  (solid 
line).  The  uniform  Peclet  number  model  is  adopted,  and  the  computations 
are  performed  using  the  steady  parabolized  approximation  on  a  stretched 
grid  with  a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the  radial 
direction  and  Nz  =  8001  points  in  the  streamwise  direction. 
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Figure  5.10:  Steady  state  temperature  profiles  at  different  streamwise  loca¬ 
tions  for  high-speed  injection  of  LP  1846  with  U  =  300m/ s,  and  two  differ¬ 
ent  inlet  temperature  To  =  278‘’J^(dashed  line),  and  To  -  298° K  (solid  line). 
The  uniform  Peclet  number  model  is  adopted,  and  the  computations  are  per¬ 
formed  using  the  steady  parabolized  approximation  on  a  stretched  grid  with 
a  stretching  parameter  a  =  6.5,  Nr  =  401  points  in  the  radial  direction  and 
Nz  =  8001  points  in  the  streamwise  direction. 


Chapter  6 

Summary  and  Conclusions 


Shear-induced  heating  of  a  liquid  monopropellant  during  high-speed  injection 
in  an  axisymmetric  nozzle  is  analyzed  numerically.  The  numerical  schemes 
are  based  on  a  finite-difference  discretization  of  the  vorticity  transport  and 
energy  equations.  Steady  and  unsteady  codes  are  applied  to  predict  peak 
temperatures  of  LP1846  during  high-speed  short-duration  injection  in  a  noz¬ 
zle  having  4mm  diameter  and  2mm  length.  When  adiabatic  wall  condi¬ 
tions  are  assumed,  computed  results  reveal  a  quadratic  dependence  of  the 
peak  temperature  on  injection  velocity.  Significant  temperature  increase,  of 
the  order  of  \0Q°K  or  more,  is  predicted  for  injection  velocities  higher  than 
200m/s.  Thus,  for  such  injection  conditions,  mixture  preignition  is  likely  to 
occur. 

Wall  heat  transfer,  modeled  in  terms  of  a  wall  heat  transfer  coefficient, 
h,  is  also  analyzed.  Three  values  h  =  20,  50,  and  SOOW/m^.^A',  which  are 
representative  of  free  convection  conditions  in  air,  forced  air  cooling  at  low 
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speed,  and  forced  liquid  cooling  at  moderate  speed,  are  selected  in  the  analy¬ 
sis.  When  the  injection  speed  is  very  high  {U  >  200m/s),  wall  heat  transfer 
dies  not  significantly  reduce  peak  temperature  predictions  even  for  high  heat 
transfer  coefficient.  In  addition,  for  low  heat  transfer  coefficient,  peak  tem¬ 
peratures  are  close  to  those  obtained  assuming  adiabatic  wall  conditions. 
Thus,  free  convection  heat  transfer  does  not  significantly  reduce  the  likeli¬ 
hood  of  mixture  ignition.  On  the  other  hand,  large  values  of  the  heat  transfer 
coefficient  may  appreciably  reduce  peak  temperatures  in  thin  walled  nozzles 
whenever  the  injection  velocity  is  not  extremely  high.  In  these  situations, 
forced  cooling  techniques  may  be  especially  tailored  in  order  to  effectively 
minimize  the  risk  of  mixture  ignition. 

The  impact  of  a  temperature-dependent  viscosity  is  also  examined.  While 
all  trends  established  using  a  constant  property  model  are  once  again  ob¬ 
served,  the  dependence  of  viscosity  on  temperature  may  significantly  affect 
shear  heating  predictions.  Computed  results  show  that,  for  moderate  injec¬ 
tion  speeds,  temperature  predictions  of  both  constant  and  variable  viscosity 
models  are  very  close.  However,  for  high  injection  speeds,  peak  temperatures 
obtained  using  a  variable-viscosity  model  are  smaller  than  those  obtained 
using  constant  property  simulation.  This  effect  is  related  to  the  decrease  of 
viscosity  and  Prandtl  number  with  increasing  temperature.  Variation  of  the 
mixture  inlet  temperature  also  affects  temperature  predictions.  In  particu¬ 
lar,  it  is  found  that  shear  heating  effects  are  more  pronounced  for  smaller 


inlet  temperature,  i.e.  when  the  inlet  Prandtl  number  is  larger. 

Finally,  it  is  worthwhile  summarizing  the  modeling  approaches  used  in 
the  present  analysis  and  the  scope  of  the  present  computations. 


(1)  In  all  multi-dimensional  simulations  discussed  here,  we  have  focused 
on  laminar  heating  in  a  small-diameter  high-speed  nozzle.  A  bound¬ 
ary  layer  of  vanishingly  small  thickness  is  simulated  by  systematically 
decreasing  the  inlet  boundary  layer  thickness  until  the  results  become 
effectively  indensitive  to  this  parameter. 

(2)  In  all  cases  considered,  a  small  constant-diameter  axisymmetric  is  as¬ 
sumed.  The  nozzle  diameter  and  length  are  kept  fixed,  d  =  4mm, 
and  R  =  2mm,  respectively.  Due  to  the  short  length  of  the  nozzle, 
the  boundary  layer  remains  much  thinner  than  the  nozzle  radius  for 
all  flow  conditions  analyzed  above.  Consequently,  no  significant  spatial 
acceleration  of  the  mixture  within  the  potential  core  was  observed.  The 
potential  contribution  of  this  effect  to  viscous  heating  of  the  mixture 
should  be  carefully  analyzed,  especially  if  nozzles  having  much  smaller 
diameter  than  those  considered  here  are  selected. 

(3)  As  mentioned  above,  the  present  analysis  of  wall  heat  transfer  is  based  on 

selection  of  a  heat  transfer  coefficient,  h.  Inherently,  the  corresponding 
models  assumes  that  the  nozzle  walls  are  very  thin;  therefore,  their 


heat  storage  capacity  is  ignored  as  is  their  contribution  to  the  overall 
thermal  resistance  to  heat  flux.  For  thick-walled  nozzles,  or  for  nozzles 
embedded  within  large  systems,  the  present  modeling  approach  should 
be  appropriately  altered. 


Part  II 

Direct  Numerical  Simulation  of 
Turbulent  Flow 
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Chapter  7 

ID  Solution  of  Turbulent  Flow 
between  Parallel  Plates 


As  discussed  in  the  Introduction,  boundary  layer  transition  is  expected  to 
occur  for  long-duration  injection  or  when  the  nozzle  length  is  large.  In  order 
to  study  shear  heating  mechanisms  within  a  transitional  or  turbulent  flow 
environment,  one  must  first  characterize  the  mean  flow  and  the  fluctuations 
around  the  mean  flow.  The  objective  of  this  chapter  is  to  construct  solutions 
for  the  mean  temperature  field,  using  well-established  empirical  correlations 
for  the  mean  velocity  field.  Temperature  fluctutations  around  the  mean  will 
be  the  focus  of  subsequent  chapters. 


7.1  The  Mean  Temperature  Distribution 

The  starting  point  for  determining  the  mean  temperature  distribution  con¬ 
sists  of  the  Reynolds-averaged  equations  of  motion.  We  focus  our  attention  on 
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fully-developed  turbulent  flow  in  a  two-dimensional  channel,  and  assume  that 
the  channel  plates  have  infinite  extent  in  the  streamwise  x-  direction,  and 
the  cross-stream  z—  direction.  For  steady,  incompressible,  constant  property 
boundary  flow,  x-momentum  and  energy  conservation  equations  are: 


du  du  1  dp  d^u 

dx  dy  p  dx  dy^ 

(7.1) 

dT  dT 

pCp{u—  +  V—)  =  kV'^T  -1-  $ 
dx  dy 

(7.2) 

where,  u  and  v  are  the  velocity  components  in  the  x  and  y  directions,  respec¬ 
tively,  p  is  pressure,  p  is  density,  Cp  is  the  specific  heat  at  constant  pressure, 
k  is  the  thermal  conductivity,  T  is  temperature,  and  $  is  the  viscous  dissi¬ 
pation  function. 


By  taking  the  time  average  of  the  energy  conservation  equation,  we  obtain 
the  following  governing  equation  for  the  mean  temperature  profile  in  a  fully- 
developed  turbulent  channel  flow: 


{_dT  _dT 


^  m 

dy 


(7.3) 


Here,  r*  denotes  the  total  shear  stress  while  qt  represents  the  total  heat  flux. 
qt  and  r*  are  respectively  given  by: 
du 


n  = 


dy 


—  pu'v' 


dT 
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and  consist  of  contributions  due  to  molecular  diffusion  and  turbulent  mixing. 

Next,  we  invoke  the  Boussinesq  analogy  for  eddy  viscosity  and  eddy  con¬ 
ductivity.  We  thus  have: 

dT 

Qt  =  —{k  +  kt)-^ 

dy 

Here,  neither  Ht  and  kf  is  a  property  of  a  fluid;  their  ratio  defines  a  dimen¬ 
sionless  group,  the  turbulent  Prandtl  number: 

Pr,  =  ^ 

Note  that  Prt  is  of  order  unity.  In  all  our  calculations,  Prt  is  taken  to  be 
0.9. 
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k  kf  —  A;(l  +  €jj/oi) 


where  cm  is  the  eddy  diffusivity  for  momentum  transfer,  while  is  the  eddy 
diffusivity  for  heat  transfer.  They  are  respectively  defined  by  the  following 
expressions: 

du  - 

P^M-^  =  -pu'v' 

dT  _ 

—  =  -v'V 
dy 

Substituting  the  above  relationships  into  the  time-averaged  energy  equa¬ 
tion,  we  get: 


d  (,  ^  .dT\  II  .  ^  ^  (du'' 

The  temperature  boundary  conditions  are: 


I 


^1  -0 
^  \v=yc  ^ 


(7.6) 


(7.7) 


Integrating  the  energy  equation  once  with  respect  to  y,  we  have 

^  ly + (I)  '‘y  (^-s) 

where  Qyj  is  the  heat  flux  at  the  wall,  which  is  determined  from  a  global  energy 
conservation  constraint.  Integrating  the  energy  equation  from  the  wall  to  the 
center  line,  we  get: 


Qw 


fVc 


dy 


(7.9) 
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Integrating  the  energy  equation  twice  with  respect  to  y,  we  have 

T-T^=  -  SlL  r 

pCp  Jo  a  +  ch 

The  above  expression  is  normalized  by  introducing  friction  velocity 


ryi  (  \  ^ 

(1  +  e„/i/)  f  —  j  dy^  1  dy,  (7.10) 


U*  =  P 


where  is  the  wall  shear  stress,  and  defining  the  following  dimensionless 


groups, 


_  {Tw  -  T)u^ 

Qw/  (pCp) 

Using  this  normalization  convention,  the  non-dimensional  temperature  dis¬ 
tribution  is  given  by: 


where. 


T+  =  - - l——dy+  +  ^ 

e, 

qw  fy-  .  fdu+V  ,  , 

^^=-1  +  ‘‘y*  (7.12) 


In  Oder  to  solve  the  above  equation,  we  need  to  estimate  the  value  of 
cm/i^,  which  depends  on  the  mean  velocity  profile.  For  simplicity,  a  two- 
layer  model  for  the  mean  velocity  distribution  is  used.  In  this  model,  the 
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turbulent  boundary  layer  consists  of  two  distinct  regions:  a  viscous  sublayer, 
in  which  u  >  cjv/,  and  a  fully  turbulent  region,  where  cm  >  i'-  In  the  viscous 
sublayer,  the  mean  velocity  is  given  by: 


(7.13) 


while  in  the  outer  layer  the  velocity  profile  follows  (Kays  and  Crawford  1980): 


Note  that  the  outer  layer  approaches 


(7.14) 


«■*■  =  2.51ny‘'’  +  5.5  (7.15) 

near  the  wall,  and  it  has  zero  slope  at  the  center  line,  i.e.  at  y  =  /i.  The  crit¬ 
ical  value  of  y"*"  for  the  two-layer  model  is  taken  to  be  12  in  our  calculations. 


Note  that  the  expression  of  cm/ in  the  sublayer  is  based  on  the  mixing  length 
model.  In  the  region  toward  the  centerline  of  the  channel  the  assumption  of 
a  constant  mixing  length,  which  is  recommended  for  the  external  boundary 
layers,  is  no  longer  appropriate.  Thus,  we  rely  on  the  empirical  equation 
proposed  by  Riechardt,  which  is  valid  for  the  entire  region  outside  of  the 
viscous  sublayer. 

Here,  h  denotes  the  distance  between  the  plate  and  central  line.  Using  the 
empirical  correlations  u'^  and  cm,  the  temperature  profile  is  determined  by 
numerically  evaluating  the  integral  on  the  right-hand  side  of  Eq.  (7.11). 


7.2  Results  and  Discussion 

The  mean  temperature  profile  for  a  fully-developed  flow  in  2D  channel  of  half¬ 
depth  h  =  2mm,  with  friction  velocity  u*  =  11  m/s  and  kinematic  viscosity 
p  =  5  X  10~^m^/s  are  obtained  using  the  approach  of  the  previous  section. 
The  corresponding  mean  velocity  u  at  the  center  line  is  about  300m/s,  and 
ranges  from  0  to  4400.  The  results  are  depicted  in  Figs.  (7.1-7.4). 

Fig.  7.1  shows  the  two-layer  model  mean  velocity  profile  in  the  range  of 
0  <  y'*’  <  4400.  Note  that  the  velocity  gradient  at  the  center  line  (y+  =  4400) 
is  zero,  as  required  by  the  symmetry  of  the  velocity  profile  across  the  chan¬ 
nel.  In  Fig.  7.2,  the  comparison  of  turbulent  and  laminar  parabolic  profiles 


for  the  same  mean  velocity  is  made.  The  velocities  are  normalized  with  the 
maximum  velocity  at  the  center  line. 

The  temperature  variations  for  different  Prandtl  numbers  are  plotted  in 
Fig  7.3.  Note  that,  the  discontinuity  of  gj  at  =  12  is  due  to  the  dis¬ 
continuity  of  at  y"''  =  12,  and  the  fact  that  the  estimates  of  e/nu  is  also 

dependent  Since  the  constant  wall  temperature  boundary  condition 

is  used  in  the  calculation,  the  wall  tempearuture  keeps  unchanged,  and  the 
highest  temperature  occurs  near  to  the  wall  due  to  the  higher  heat  generation 
there. 

Finally,  the  computed  temperature  profile  is  compared  in  Fig.  7.4  with 
that  of  laminar  entrance  flow.  The  plots  indicated  that  a  higher  peak  tem¬ 
perature  is  predicted  for  fully-developed  turbulent  flow,  and  this  effect  is 
attributed  to  the  fact  that,  in  the  turbulent  flow  region,  the  viscous  sub¬ 
layer  is  spatially  uniform  and  has  a  thickness  which  is  significantly  smaller 
than  that  of  the  laminar  boundary  layer  in  the  entrance  region. 


Figure  7.1:  Turbulent  mean  velocity  distribution  in  wall  coordinates.  The 
calculation  is  performed  with  maximum  mean  velocity  Umax  =  300m/s,  half 
channel  depth  h  =  2mm,  viscosity  u  =  5.0e“®. 
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Figure  7.2:  Comparison  of  turbulent  and  laminar  parabolic  profiles  for  the 
same  mean  velocity.  The  velocities  are  normalized  with  the  maximum  tur¬ 
bulent  velocity  at  the  center  line.  The  turbulent  velocity  profile  is  computed 
with  maximum  mean  velocity  u^ax  =  300m/s,  half  channel  depth  h  =  2mm, 
viscosity  u  =  5.e“®. 


Figure  7.3:  Turbulent  mean  temperature  profile  with  four  different  Prandtl 
numbers,  7,  25,  64  and  100  in  wall  coordinates  and  constant  wall  temperature 
condition.  The  calculations  are  performed  with  maximum  mean  velocity 
'^max  “  300m/ s,  half  channel  depth  h  =  2mm,  and  viscosity  u  =  5.e“®. 


Figure  7.4:  Dimensional  turbulent  mean  temperature  profile  (maximum 
mean  velocity  u^ox  =  300m/s)  and  laminar  entrance  temperature  profile 
(maximum  velocity  is  300m/s)  with  Pr  =  109,  half  channel  depth  h  =  2mm, 
and  viscosity  v  =  5.e“®.  Also,  constant  wall  temperature  condition  is  im¬ 
posed  at  both  cases. 


Chapter  8 

Homogeneous  Isotropic 
Turbulence 

8.1  Introduction 

As  discussed  earlier,  characterization  of  the  effects  of  shear-induced  heating 
in  transitional/turbulent  flows  necessitates  a  study  of  unsteady  temperature 
fluctuations  around  their  mean  values.  In  the  present  chapter,  estimates 
are  obtained  for  the  simplified  case  of  a  statistically-steady  homogeneous 
isotropic  turbulent  field.  The  primary  motivation  behind  the  present  exer¬ 
cise  is  to  derive  approximate  estimates  for  the  amplitude  of  spatial  tempera¬ 
ture  fluctuations,  and  to  gain  insight  into  the  phenomena  which  govern  their 
behavior.  To  this  end,  we  shall  select  high  values  for  the  kinetic  energy  dis¬ 
sipation  rate  which  are  characteristic  of  high-speed  flows  in  thin  channels. 
The  evolution  of  the  isotropic  velocity  and  temperature  distributions  will  be 
computed  using  direct  numerical  simulation  of  the  vorticity  transport  and 
energy  equations.  The  results  will  then  be  contrasted  in  the  following  chap>- 
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ter  with  simulations  of  turbulent  channel  flow. 


8.2  Formulation 


We  assume  an  incompressible,  constant-density  fluid  with  constant  proper¬ 
ties.  The  motion  of  the  fluid  is  governed  by  the  conservation  equations  for 
mass,  momentum  and  energy.  In  a  vorticity-based  formulation,  the  governing 
equations  are  expressed  as: 


0  UJ  ^  o  ^ 

-b  V  X  (a;  X  u)  =  u)  -fVx  / 

u=  -Vx  u 
DT 

pCp—  =  kV^T  +  n^ 


(8.1) 

(8.2) 

(8.3) 


where  u  is  the  velocity  vector,  u  is  the  vorticity,  T  is  temperature,  p  is  den- 

— y 

sity,  t  is  time,  f  is  the  force  acting  on  the  fluid,  Cp  is  speciflc  heat  at  constant 
pressure,  k  is  the  thermal  conductivity,  u  is  the  kinematic  viscosity,  p,  is  the 
dynamic  viscosity,  ^  w  •  V  is  the  material  derivative,  and  $  is  viscous 

dissipation  function.  The  above  vorticity-velocity  formulation  has  become, 
to  a  growing  number  of  people,  an  attractive  alternative  to  formualtions  in 
primitive  variables.  This  interest  (Guevremont  et  al.  1990;  Gresho  1991)  is 
mainly  linked  to  easier  treatment  of  boundary  conditions  since  the  pressure 
is  no  longer  part  of  the  solution. 


Let  b=  Vx  /,  /=  (/i,/2,/3),  h=  (&i,62,^>3)  in  a  Cartesian  coordinate 
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system  {x,y,z)  —  {xi,X2,xz).  In  component  form,  the  vorticity  forcing  term 

is  given  by: 

.  ^  dfz  df2 

^  dy  dz 

^  dz  dx 

^  dx  dy 

and  the  governing  equations  are  expressed  as: 


du>i  dS3  dS2  _2  ,  , , 

(8.4) 

duj2  dsi  ds3 

(8.6) 

du}3  ds2  dsi  2  ,  .0 

dt'^  dx  =  +  « 

(8.6) 

2  _  9^2  du)3 

dz  dy 

(8.7) 

ax  dz 

(8.8) 

dy  dx 

(8.9) 

dy  dz  c„ 

(8.10) 

Here, 


Si  =  U3U2  -  U2OJZ 
$2  =  U1U3  —  U3U1 
S3  =  U2a;i  —  U1U2 

and  the  viscous  dissipation  function  $  is  given  by: 


(8.11) 
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+ 


'dui  du2 
dy  dx 


+ 


fdui  duz 
\dz  dx 


+ 


'  du2  duz 
.dz  dy  ^ 


For  a  statistically-steady,  isotropic  turbulent  flow,  we  further  assume  that 
periodic  boundary  conditions  hold  along  the  coordinate  directions  for  the 
velocity,  vorticity,  and  temperature  fields.  Thus,  these  variables  admit  a 
Fourier  representation  in  terms  of  the  coordinate  variables  x,  y  and  In 
Fourier  space,  the  governing  equations  take  the  following  form: 


+  k  P)wi  =  ik2S2  -  ikySz  +  h 

,  d  I  .9.  ^ 

(^  +  ^1  A:  I  )^2  =  -  ikzSi  +  62 

^  ^ 

(^  +  1^1  k  P)a'3  =  ikySi  —  ikxS2  +  bz 

I  /c  pui  =  ikzU}2  —  ikyU)z 
I  k  pU2  =  ikxi^z  —  ikzCji 
I  k  I  U3  —  zkyCOi  ikx(-02 


d 


(—  +  a\  k  P)T  +  ikxUiT  +  ikyU2T  +  ik^uzT 


^dt  »  -  cp 

— f 

where,  k=  {kx,  ky,  k^)  is  the  wavenumber  vector,  and 


k 


=  ^ki+ki+ki 


(8.12) 

(8.13) 

(8.14) 

(8.15) 

(8.16) 

(8.17) 

(8.18) 

(8.19) 


is  the  modulus  of  I  k 


8.3  Numerical  Schemes 

Flowfield  simulation  is  performed  using  a  pseudo-spectral  formulation  of  the 
governing  equations.  Briefly,  the  evolution  of  the  Fourier  coefficients  for 
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vorticity  and  temperature  are  computed  by  approximating  non-linear  source 
terms  using  spectral  collocation  derivatives  and  numerically  integrating  the 
resulting  equations. 

In  the  numerical  implementation,  we  start  with  the  vorticity  transport 
equation.  We  rely  on  factorization  of  the  viscous  diffusion  term,  and  recast 
the  equations  of  motion  of  the  Fourier  coefficients  of  vorticity  as: 


^{ikzS2  -  ikySs  -1-  ^1) 

(8.20) 

+  (2) 

(8.21) 

-  ihsi  +  h) 

(8.22) 

The  Fourier  coefficients  Si, 

S2,  and  S3  are  approximated  using  spectral  collo- 

cation  derivatives  and  the  discrete  system  is  integrated  using  the  second-order 

Adams-Bashforth  scheme. 

The  resulting  discrete  evolution  equations 

are: 

g-n+\  _  g-,/|tpAt^n 

J=0 

At 

-  iky^^-^)  -f 

(8.23) 

^n+1  _  g-HfcpAf^n 

j=0 

At 

-  ikzsl'^~^)  +  e“''l 

(8.24) 

_  g-./|fepAt— n 

i=o 

At 

(8.25) 

In  equations  (8.23-8.25),  i  =  v^,  Po  =  3/2,  A  =  -1/2,  and  At  is  time 
step. 

Recall  that  the  Fourier  transform  of  f{x)  is  defined  by: 

1  r+oo 

^xlfix)]  =  F(0  =  fix)e~^^^dx  (8.26) 

and  that  the  inverse  Fourier  transform  of  F{^)  is  given  by: 

n  1  r+OO 

fix)  =  F-‘[F(?)1  =  (8.27) 

where  Fx  and  Fx  ^  denote  the  forward  and  inverse  Fourier  transforms,  re¬ 
spectively.  By  considering  terms  of  the  form  F~^[e~°'^^ Fx[f  {x)]],  where  a 
is  a  positive  number,  it  is  easy  to  verify  using  the  definitions  of  Fx  and 
Fx^  that  if  the  function  f{x)  is  real- valued  then  this  is  also  the  case  for 
^F^[^~°‘^^Fx[f{x)]].  A  similar  observation  also  holds  for  terms  of  the  form 
Fx-^[iCe-<^^Fx[f{x)]]. 

In  the  computations,  we  take  advantage  of  the  above  observations,  and 
also  exploit  the  linearity  of  Fourier  operators  and  the  fact  that  transforms 
and  inverse  transforms  along  different  coordinate  directions  commute.  By 
doing  so,  we  arrive  at  the  following  discrete  system  of  evolution  equations: 

J=0 

-  {F-^e-‘'>‘"^^Fx){F;^e-<^HkyFy){FF'^e-’'^^^^^ 
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j=0 

+  At(F-h-'''‘"^^F,){F;^e-‘'’^l^*Fy){F-^e-^^^^^^  (8.29) 

0;?+'  =  {F-h-‘''‘"^^F,){Fy-^e-'''‘"y^^Fy){F-^ 

j=0 

[  {F-^e-^’^'^^F^){Fy^e-'''‘y^HkyFy){F-h-''^^^^^ 

-  {F-h-^'‘"^Hk^F^){Fy-h-‘''^y^^Fy){F-h-^^^^^ 

+  Ai(F-'e-'''='^‘F^)(F-ie-''*^vA‘Fj,)(F7'e-'''='^‘F^  (8.30) 

The  advantage  of  the  above  system  is  that  it  yields  the  discrete  evolution 

of  the  vorticity  field  in  physical  space  using  a  spectral  collocation  scheme 
which  involves  real  Fourier  transforms  only.  This  simplifies  the  numerical 
implementation  and  results  in  significant  computational  savings  over  com¬ 
plex  FFTs. 

Once  the  vorticity  field  is  updated,  Eq.  (8.2)  is  inverted  in  order  to  de¬ 
termine  the  velocity  field.  In  the  computations,  we  first  consider  the  velocity 
components  ui  and  U2,  as  well  as  the  combination, 

Uc  =  ui  +  iu2  (8.31) 

Based  on  the  updated  values  of  Wi,  u)2,  ^3,  we  form  the  source  terms 

_  dLiJ2  du}^ 

dz  dy 
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_  5^3  diVi 
~  dx  dz 

using  spectral  collocation  derivatives,  and  define  the  combination, 


rc  =  ri  +  1X2 


(8.32) 


Using  Eqs.  (8.15)  and  (8.16),  the  Fourier  transform  of  Uc  is  determined  from: 


0 


(8.33) 


with  Uc  —  0  for  k=  0.  The  components  ui  and  U2  are  then  solved  together 
by  inverting  Eq.  (8.33) ,  at  the  cost  of  one  complex  FFT. 


Once  ui  and  U2  are  determined,  the  velocity  component  U3  is  obtained 

by  using  continuity  equation  and  the  definition  of  vorticity,  Vx  u.  In 

Cartesian  coordinates,  we  have: 

du3  du2 
dy 


Ui 


U)2  = 


dui 

dz 

dU2 

dx 


dz 

du3 

dx 

dui 

dy 


(8.34) 

(8.35) 

(8.36) 


From  Ml  and  u)2,  we  find  using  Eq.  (8.35),  and  its  Fourier  transform 
with  respect  to  x.  Suppose  that  ^  =  Ur  +  ioj,  and  m5  =  uir  +  then 
ikxU3  =  Or  +  ihi.  When  fca,  ^  0,  M3  can  be  obtained  through  the  spectrum 
of  Assuming  M3  vanishes  at  (0,0,0),  we  have  from  the  definition  of  the 
Fourier  transform: 


NI2-1 

M3(0,0,0)=  Y.  U3{kxM 

kx=-Nl2 


(8.37) 


Note  that  for  1-D  real  FFT,  the  real  part  of  the  Fourier  coefficients  is  sym¬ 
metric  with  respect  to  wave  number,  while  the  imaginary  part  of  the  coeffi¬ 
cients  is  antisymmetric  with  respect  to  wave  number,  i.e.  -Re(w3(A:i,0,0))  = 

■Re(w3(— 0}  0)),  Im{u3{kx,0,0))  =  — W3(— Aix,  0, 0).  Thus,  the  above  expres¬ 
sion  reduces  to 

N/2-1 

ui(0,0,0)  =  -2  Re(u^(kx,0,0))  (8.38) 

kx~l 

and  one  inverse  FFT  will  yield  U3(a;,  0, 0),  i.e  velocity  component  along 
the  ar-axis. 


Next,  from  M3(a;,  0, 0),  we  can  solve  for  uz{x,y,0),  i.e  the  value  of  in 
the  x-y  plane.  We  now  form  using  Eq.  (8.34),  and  its  Fourier  transform 
with  respect  to  y  to  obtain  an  expression  of  the  form,  ^  =  hr  +  ih-  Then, 
'^3  =  bi/ky-\-ibj./ky  if  ky  ^  0.  The  coefficient  of  the  zeroth  mode  is  determined 
using  the  definition  of  the  Fourier  transform;  we  have: 

N/2-1 

M3(a;,0,0)=  ^  ui{x,ky,0)  (8.39) 

ky  =  -N/2 

Taking  advantage  of  the  conjugate  properties  of  the  Fourier  coefficients,  we 
obtain; 

N/2-1 

U3{x,ky  =  0,0)  =  u3{x,y  =  0,z  =  0)-2  Y,  U3{^,ky,0)  (8.40) 

ky=l 

where  U3(a:,0,0)  is  known  from  the  previous  step.  Then,  the  inverse  FFT 
will  give  the  value  of  in  x-y  plane. 


Ill 


Finally,  we  consider  the  continuity  equation  to  compute 
dus  dui  du2. 

and  apply  a  similar  procedure  to  that  of  the  previous  two  steps  in  order  to 
determine  the  Fourier  coefficients  ui{x,y,k^).  The  complete  distribution  of 
U3  is  then  determined  using  ID  inverse  Fourier  transforms  for  points  in  the 
x-y  plane. 


The  new  temperature  field  is  obtained  in  the  last  step  of  the  numerical 
scheme.  The  energy  equation  is  treated  in  a  similar  fashion  to  the  integra¬ 
tion  of  the  vorticity  transport  equation.  Specifically,  we  rely  on  spectral 
collocation  derivative  for  the  convective  terms  and  viscous  dissipation  func¬ 
tion  and  on  factorization  of  diffusive  terms.  In  conjunction  with  this  spatial 
discretization  scheme,  we  use  the  second-order  Adams-Bashforth  scheme  for 
the  convective  terms,  and  treat  the  dissipation  term  in  an  implicit  fashion 
using  a  Crank-Nicolson  formulation.  The  resulting  discrete  temperature  evo¬ 
lution  equation  is: 


1^— ] 


+4-.+1J  _ 

j=0 


(8.42) 


where  P  =  V.{u  T). 


8.4  Forcing 


In  order  to  maintain  a  statistically  steady  field,  a  forcing  function  is  used  on 
the  right-hand  side  of  the  vorticity  transport  equation  with  forcing  applied  at 
lower  modes.  Specifically,  all  Fourier  modes  with  wavenumber  components 
equal  to  0  or  1  are  forced  with  a  constant  amplitude  a/,  but  with  a  random 
phase  6.  Both  the  amplitude  and  phase  are  independent  of  t.  Note  that  the 
components  of  k  take  only  integer  values,  since  the  space  period  is  27r.  Such 
a  force  field  can  be  synthesized  as 

111 
fj{x,y,z)=  ^  ^  ^ 

kx^O  Ajy~0  kz^O 

and  the  three  components  of  force  field  are  generated  independently.  As  will 
be  shown  later,  the  kinetic  energy  dissipation  rate  does  not  achieve  constant 
steady-state  values  using  the  present  forcing  approach.  Thus,  a  statistically- 
steady  state  will  be  reached  in  the  sense  that  the  dissipation  rate  undergoes 
small-amplitude  fiuctuation  around  a  mean,  constant  value. 

To  reach  a  statistically  steady  state,  energy  loss  through  viscous  dissipa¬ 
tion  had  to  be  balanced  with  the  energy  injection  through  large-scale  foring 
function.  Since  we  wish  to  fix  the  value  of  the  kinetic  energy  dissipation  rate, 
€,  the  magnitude  of  the  forcing  field  is  adjusted  during  the  calculations.  In 
particular,  the  coeflScients  of  the  forcing  field  are  adjusted  every  several  time 
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steps  so  that: 

^jrudy  =  t 

where  V  is  the  volume  of  the  compuational  domain. 


(8.44) 


8.5  Initial  Conditions 

A  random  velocity  field  u  with  —  |  spectrum  is  synthesized  by  superposing 
Fourier  modes  with  random  phases  6,  uniformly  distributed  in  the  range  of 
[0,27r], 

Ui{x,  y,z)=  {kl  +  kl  +  (8.45) 

hx jhz 

The  three  velocity  components,  Uj  z  =  1, 2, 3,  are  generated  independently. 
Since  the  resulting  velocity  field  may  not  be  divergence  free,  a  correction 
(or  projection)  step  is  necessary  to  ensure  suitable  initial  conditions.  There 
exist  several  approaches  which  can  be  performed  on  the  synthesized  field  Ui 
to  render  it  divergence  free.  The  method  we  have  used  is  based  on  finding 
the  vorticity  of  the  synthesized  field,  then  reconstructing  the  velocity  field 
by  inverting  Eq.  (8.2).  This  guarantees  that  the  reconstructed  field  satisfy 
the  divergence-free  condition. 


As  mentioned  earlier,  we  wish  to  prescribe  energy  dissipation  rates  that 
are  characteristic  of  high-speed  flows  in  pipes  or  channels  with  small  cross- 
sectional  area.  As  a  representative  example,  we  consider  a  pipe  of  diameter 


D,  with  mean  flow  velocity  Um-  Based  on  these  ‘input’  values,  we  determine 
the  Reynolds  number, 


(8.46) 


where  u  is  the  kinematic  viscosity.  Next,  we  take  advantage  of  well-known 
correlations  for  the  friction  factor  (Incropera  &  DeWitt  1990), 

f  /  =  0.316Rec^/'‘  ReD<2x  10^ 

\  f  =  0.184Re^^/^  ReD>2x  10^ 

use  the  definition  of  the  friction  factor. 


/ 


—{dp/dx)D 

P</2 


(8.48) 


and  the  force  balance  across  the  area  of  the  pipe  to  obtain  the  friction  veloc¬ 
ity: 


D1  dp 
A  pdx 


(8.49) 


For  turbulent  pipe  flow,  the  viscous  dissipation  rate  c  can  be  approximated 
as 


er^ul/Ky  (8.50) 

where  y  is  the  distance  from  the  wall  and  k  is  the  von-Karman  constant.  In 
the  calculations,  we  shall  choose  values  of  e  corresponding  to  y  =  12u fu*,  i.e. 
the  viscous  dissipation  rate  at  the  viscous  sublayer. 


8.6  Results  and  Discussion 


We  present  results  of  three  simulations,  performed  with  different  resolution 
grids.  The  first  simulation  is  performed  on  a  grid  with  x  Ny  x  = 
64  X  64  X  64  mesh  points,  while  the  other  two  are  computed  with  a  grid  hav¬ 
ing  Nx  X  Ny  X  Nz  =  96  X  9Q  X  96  points.  For  all  three  cases,  the  prescribed 
viscous  dissipation  rate  e  =  1.25  x  lO^m^/s^.  Following  the  discussion  of  the 
previous  section,  the  chosen  value  of  e  corresponds  to  a  high-speed  flow  with 
mean  velocity  Um  =  300m/s,  viscosity  u  =  5.0  x  lQ~^m^/s,  in  a  pipe  of  inter¬ 
nal  diameter  D  =  4mm.  Based  on  the  prescibed  value  of  e,  the  Kolmogorov 
microscale  rj  =  ~  0.562/xm. 

The  conditions  of  the  simulations  are  listed  in  Table  8.1  which  in  par¬ 
ticular  provides  values  of  the  Kolmogorov  microscale,  77,  the  integral  scale, 
Lc,  the  Prandtl  number,  Pr,  as  well  as  the  minumum  allowable  values  of 
the  normalized  Kolmogorov  scale,  r]min,  and  the  normalized  Batchelor  scale 
^min-  fa  estimating  these  quantities,  which  are  governed  by  the  resolution 
of  the  grid,  we  have  related  the  Batchelor  scale  to  the  Kolmogorov  scale  by 
^  =  y/Prrj,  and  the  criterion  K  x  max{fj,i)  >  1.5  for  the  flow  to  be  well 
resolved.  Here,  K  =  N/2  is  the  normalized  wavenumber  corresponding  to 
wavelength  equal  to  one  grid  size.  Note  that,  since  resolution  requirements 
increase  significantly  with  increasing  Prandtl  number,  we  have  restricted  the 
computations  to  the  range  Pr  <  13. 


Table  8.1 


Case 

Vmin 

V 

Lc 

Lch 

Pr 

1 

0.0468 

0.1146 

0.562//m 

31//m 

55.3 

2 

0.03125 

0.07654 

0.562)um 

46fj,m 

81.8 

6 

3 

0.03125 

0.1148 

0.b62fj,m 

Slum 

55.3 

13.5 

The  computations  are  carried  out  for  simulation  times  that  are  long 
enough  for  statistically  steady  conditions  to  be  reached.  To  verify  that  this 
is  in  fact  the  case,  we  plot  in  Fig.  8. 1-8.2  the  evolution  of  the  dissipation 
rate,  e,  and  the  root-mean-square  velocity,  Urms-  The  evolution  of  the  energy 
spectra  for  cases  1-3  is  depicted  in  Figs.  8.3-8.5.  The  figures  indicate  that 
the  root-mean-square  velocity  undergoes  small-amplitude  fluctuations,  and 
that  the  energy  spectra  taken  at  large  simulation  times  are  nearly  identi¬ 
cal.  Thus,  statistically-steady  conditions  are  reached  during  the  simulations, 
and  this  is  not  surprising  since  the  simulation  times  correspond  to  several 
eddy  turnover  times.  For  the  presently  selected  value  of  the  dissipation  rate, 
e  ~  1.25  X  the  eddy  turnover  time  is  about  0.38/is  for  case  1, 

0.42/is  for  case  2,  and  0.40/xs  for  case  3. 

The  primary  results  of  the  temperature  simulations  are  summarized  in 
Figs.  8.6-8.8,  which  respectively  depict  the  temperature  spectrum  at  dif¬ 
ferent  time  steps  for  the  cases  1-3.  The  figures  indicate  that,  after  a  few 
turnover  times,  the  temperature  spectrum  also  reaches  a  statistically  steady 
state.  The  amplitudes  of  the  Fourier  coeflicients  remain  quite  small,  of  the 


order  of  1°K  or  less.  This  leads  us  to  expect  that  spatial  temperature  fluctu¬ 
ations  also  remain  small.  In  order  to  verify  this  expected  trend,  we  compare 
in  Fig.  8.9  the  evolution  of  the  peak  within  the  domain  to  that  of  the  average 
temperature.  We  start  from  a  uniform-temperature  initial  condition,  so  that 
the  peak  and  mean  temperatures  coincide  at  time  t  =  0.  Due  to  the  (nearly) 
constant  dissipation  rate,  the  mean  temperature  rises  linearly,  at  a  rate  of 
about  1.6^K/fis.  Despite  the  high-rate  of  energy  injection,  and  the  rapid 
increase  of  the  mean  temperature  of  the  fluid,  the  peak  temperature  remains 
within  less  than  two  degrees  from  the  average.  This  trend  is  consistent  with 
the  results  for  the  temperature  spectrum  which  exhibits  amplitude  levels  of 
the  same  order  approximately. 

We  also  note  that  the  amplitude  of  the  spatial  fluctuations  in  tempera¬ 
ture  are  weakly  dependent  on  the  Prandtl  number.  This  observation,  may 
be  verified  by  comparing  the  range  of  amplitudes  in  the  temperature  spec¬ 
tra  in  Figs.  8. 6-8.8,  should  be  contrasted  with  the  quasi-one-dimentional 
predictions  in  the  previous  chapter  which  refiected  a  linear  dependence  on 
Pr.  Thus,  although  the  present  simulations  are  restricted  to  low  Reynolds 
numbers,  the  present  results  support  the  notion  that  in  isotropic  turbulence 
the  regions  of  high  energy  dissipation  also  experience  vigorous  mixing,  and 
that  the  turbulent  mixing  of  the  fluid  substantially  reduces  otherwise  sharp 
thermal  gradients.  In  the  following  chapter,  this  intuitive  picture  will  be 
re-examined  in  light  of  simulations  of  turbulent  channel  flow. 


lO’i 


run5-a 


spectrum  for  case 


Chapter  9 

Turbulent  Channel  Flow 


9.1  Introduction 

In  this  chapter  we  examine  the  effects  of  shear-induced  heating  in  a  plane 
channel  flow.  As  indicated  below,  we  adapt  the  vorticity-based  formulation 
of  the  previous  section  to  periodic  channel,  and  rely  on  direct  numerical  sim¬ 
ulation  of  the  governing  equations  to  characterize  the  temperature  field  in  a 
transitional/ turbulent  flow  environment.  The  predictions  are  then  contrasted 
with  earlier  results  for  isotropic  turbulent  flow. 


9.2  Formulation 

The  physical  formulation  is  identical  to  that  of  the  previous  chapter.  Specif¬ 
ically,  we  assume  an  incompressible  fluid  with  constant  physical  properties, 
and  adopt  a  vorticity-based  formulation  of  mass,  momentum  and  energy  con¬ 
servation  equations. 
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In  order  to  describe  the  set-up  of  the  computations,  we  recall  from  the 
previous  chapter  the  governing  equations  in  component  form.  We  have: 


du^i  dS3  dS2  _2 

at  +  Sj, 

(9.1) 

du}2  .  dsi  ds3  2 

(9.2) 

du3  ds2  dsi  2 

dt  dx  dy  ^ 

(9.3) 

2  9u2  du)s 

(9.4) 

V  V  = 

dx  dz 

(9.5) 

.  ^  a.,,  8u., 

"  dy  dx 

(9.6) 

dt  dx  dy  dz  Cp 

(9.7) 

where  {u\,U2,uz)  denotes  the  velocity  vector,  (wi,  a;2,  0^3)  the  vorticity. 


5x  = 


S2  = 


53  =  U2U)i  —  U1U2 


t  is  time,  u  the  kinematic  viscosity,  a  the  thermal  diffusivity,  Cp  the  specific 
heat  at  constant  pressure,  and  $  the  viscous  dissipation  function. 


The  governing  equations  are  solved  in  periodic  channel  of  height  2d.  The 
channel  walls  are  assumed  to  be  flat  with  normal  along  the  2;-direction.  The 
mean  flow  is  oriented  along  the  rc-direction  and  y  is  the  spanwise  direction. 
We  assume  that  the  velocity  and  temperature  field  are  periodic  along  both 
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the  X  and  y  directions.  No  slip  boundary  conditions  are  imposed  at  solid  walls 
are  imposed,  and  heat  losses  at  the  walls  are  modeled  using  a  prescribed  heat 
transfer  coefficient.  At  the  solid  walls,  z  =  0  and  2  =  2d,  we  have: 

Ml  =  U2  =  W3  =  0 

dT 

-k—  =  h{T^  -  Teo) 

where,  d  is  the  half  depth  of  the  channel,  h  is  heat  transfer  coefficient,  n  is 
normal  to  the  wall,  is  wall  temperature,  and  Too  is  the  ambient  temper¬ 
ature  outside  channel.  In  addition  to  the  above  conditions,  the  definition  of 
the  vorticity  is  also  enforced  at  the  solid  boundaries.  As  explained  below,  the 
mean  vorticity  at  the  walls  is  related  to  the  mean  pressure  gradient,  which 
is  treated  as  a  parameter  of  the  problem. 


9.3  Normalization 


The  numerical  scheme  discussed  in  the  following  sections  solves  a  normalized 


form  of  the  governing  equations.  Variables  are  normalized  using  the  appro¬ 
priate  combination  of  the  the  fluid  density,  p,  the  channel  half  depth  d,  the 
centerline  velocity  Uc  of  the  undisturbed  flow,  and  the  ambient  temperature 
outside  of  the  channel  T^o-  For  this  choice  of  characteristic  density,  length, 
velocity  and  temperature  scales,  the  normalized  governing  equations  become: 


^  ^  ^ 

dt  dy  dz 


=  — vV 

Re  ^ 


du}2  dsi 
dt  dz 


dx  Re 


(9.8) 

(9.9) 


130 


dT  d{u{r)  diu2T) 

dt  dx  dy 

where 

Re  =  ^ 

is  the  Reynolds  number, 


a 


the  Prandtl  number,  and 


Ec  = 


CpToo 


dujj 

dt 


+ 


djuzT) 

dz 


dS2  dsi  1  2 

dx  dy  Re" 

(9.10) 

^ 

dz  dy 

(9.11) 

2  9uj3  dui 

\  Uo  = - 

dx  dz 

(9.12) 

^  dy  dx 

(9.13) 

1  1  —2m  •E'C, 

-  =  V^T  -1 - $ 

RePr  Re 

(9.14) 

(9.15) 

(9.16) 

(9.17) 

the  Eckert  number. 


9.4  Numerical  Scheme 

The  numerical  scheme  used  is  based  on  a  mixed  pseudo-spectral  finite-difference 
discretization  of  the  governing  equations.  We  use  Fourier  expansions  in  x  and 
y  directions,  and  rely  on  second-order  centered  difference  in  the  cross-stream 
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2  direction.  By  Fourier  transforming  the  governing  equations  in  both  x  and 
y,  we  obtain  the  following  governing  system  for  the  Fourier  coefficients: 


doJi  5s2 


du>2  ..  _  dsi  1  I  ,o  _ _ 

ar  "  *  I  + 

—  +  ^A:a.S2  -  zA:j,Si  = —{-\  k  \^  +  —)u}3 


Re' 


dz^‘ 


duj-i 


(-|  k  P  +  ^)“2  =  - 


5a;i 

dz 


dT  .  ^  ^ 

h  ik^Tfi  +  ikyTf)  + 


(-|  k  r  +  -  ikxu;2 

fr  =  ^-(-1  J  P  +  —)f  +  — $ 

i?ePr^  '  ^  dz'^^  ^  Re 


(9.18) 

(9.19) 

(9.20) 

(9.21) 

(9.22) 

(9.23) 

(9.24) 


where 


Ta  =  UiT 
Tb  =  U2T 
Tc  =  U3T 


The  boundary  conditions  for  the  above  system  are: 


Wi  |i=  Ui  0 

U2  |i=  ^  |n=  0 


“3  |i—  U3  |Ar=  0 
^  li.A,=  0  (A;.  =  ky  =  0) 

Wi  li,ji=  (I  A:  I  #  0) 
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I  ~az  ll.^~  Re^  —  ky  —  0) 

\  ^2  |i,Ar=  ^  (I  A:  I  7^  0) 

^  |l,7V=  0 
3T 

\l,N=  Kit  |i,;v  -1)  (9.25) 

Here,  Nu  =  hd/k  is  the  Nusselt  number,  and  N  denotes  the  number  of  fi¬ 
nite  difference  points  in  the  z  direction.  Note  that  the  equation  of  ^3  is  not 
coupled  with  other  equations,  so  that  its  time  integration  may  be  carried  out 
independently  of  the  others.  On  the  other  hand,  aJi  and  U2  are  coupled  at 
the  boundaries  with  U2  and  ui.  To  overcome  the  difiiculties  associated  with 
this  coupling,  a  boundary  Green’s  function  technique  is  used. 

Numerical  integration  of  the  vorticity  transport  and  energy  equations  is 
based  on  linear  multi-step  methods.  For  diffusion  terms,  an  implicit  treat¬ 
ment  is  used  based  on  the  second-order  Crank-Nicolson  scheme.  We  use  the 
second-order  Adams-Bashforth  scheme  for  the  convective  terms  in  the  vor¬ 
ticity  transport  equation,  and  rely  on  spectral  collacation  derivative  for  the 
convective  terms  in  the  energy  equation.  Following  Daube  (1992),  a  stag¬ 
gered  grid  is  employed  in  2;  direction.  The  grid  points  for  the  different  field 
quantities  are  defined  as  follows: 

•  ui,  U2  and  a;3  are  computed  at  nodes  {iAxJAy,  {k  +  ^)Az)  for  i,j  = 

k  =  l,---,N  -1. 

•  U3,  uji,  U2  and  T  are  computed  at  nodes  {iAx,jAy,kAz),  for  i,j,k  = 
!,•••, AT. 


133 


At  each  interior  point,  standard  second-ordered  differences  have  been  used 
for  both  first  and  second  derivatives.  For  nodes  adjacent  to  the  channel  walls, 
a  special  treatment  is  employed.  Specifically,  for  velocity  component  ui,  the 
second  derivative  is  approximated  using 
4  . 

This  approximation  comes  from  the  usual  centered  diflference  over  the  nodes 
{hji  |)>  (*5  §))  and  {i,j,  |)  incorporated  with  the  extrapolation  of  ui  at  the 

fictitious  point  {i,j,  1)  outside  of  the  computational  domain 

i  =  5  —  6u,  3  +  8ui  i) 

1,2  Q  V  1,2  ^12  ’  ' 

Similarly,  the  first  derivative  of  ui  on  the  nozzle  wall  A;  =  1  is  approximated 


using: 


.dui 
^  dz 


;(-Ui  5  -f- 9ui  I  -  8Ui,i) 


Thus,  derivatives  in  the  2-directions  are  approximated  using  second-order 
differences  at  both  the  interior  and  boundary  nodes. 


Implementation  of  the  numerical  scheme  is  summarized  as  follows.  The 
Fourier  coefficient  for  the  2-component  of  vorticity,  uj^  is  updated  by  inverting 
the  discrete  Helmholtz  equation: 

=  RWP  (9.26) 

with  Dirichlet  boundary  conditions: 

/  =  0 
1  ^3^^^  =  0 


(9.27) 


Here,  H  denotes  the  Helmholtz  operator, 


Hisr' = 


+  (i  +  ^4rl  t  P 


2Red.z^  *-5  2Re 

■  ^  r)“3r*i+  ' 


^3"^3 


2i2eAz2^  •’i+t  2i2eA^2“'^i+t 

and  i2PF  the  explicit  source  term, 

~  ~  i^xS2i^i  +  ikySn^i  -  ^^1  k  P^”+i 

1  ^"+1  “  2a5i^^i +W3”_i 
2i2e 

Next,  we  solve  for  u{  and  U2  by  inverting  the  coupled  system: 


(9.28) 


L2^"11  =  Q2W{u2'^-^^) 

*“2 


with  boundary  conditions 

/  5ir‘ = 0 

\  SIS,-"'  =  0 

155,  =  ^  if|fc|#0 
lS?Ar  =  2|f  it|*|#0 

/  ^  =  if  1*1  =  0 
I  If  =  -R<=iU  ifW  =  o 

Here,  H2  denotes  the  Helmholtz  operator. 


1 


2ReAz^ 


At  2Re' 


2ReAz^ 


rr-n+l 

^2i+i 


(9.29) 

(9.30) 

(9.31) 

(9.32) 

(9.33) 


(9.34) 
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and  R2W  the  explicit  source  term, 


R2W{cr2'i) 


^  Ru-\ 

Ai  Az 

1  l^p^n  I  1  ^"+1  -  25^"  +  ^7-1 


Meanwhile,  L2  represent  a  discrete  “Laplacian”  operator  defined  by: 


(9.35) 


(9.36) 


rT-n+l  .*-7-71+1 

Qmch'Xi)  = 


(9.37) 


is  the  corresponding  source  term. 

For  non  zero  modes,  boundary  Green’s  functions  are  employed  to  decouple 
the  equations  of  Ui  and  A2.  The  boundary  Green’s  function  technique  exploits 
the  linearity  of  the  elliptic  operators  by  first  solving  the  auxiliary  system: 


r  =  i22W(a5^n 
\  5}^”+'  =  ^^+1  =  0 

r  T2«i”4  =  g2iF(^ri) 

1  Sir'  =  =  0 

and  then  expressing  the  solution  of  the  full  system  as: 

f  ar2r'=i5;r')+ELiai«? 

\u^:i=si7:i+EUa.f^ 

Here,  gf  and  are  the  solutions  of  the  elementary  problems: 

'  H29^  =  0  {k  =  l,2) 

^  =  1  9n  =  ^ 
i5?  =  0  gl  =  l 


(9.38) 


(9.39) 


(9.40) 


(9.41) 


L2ft  =  0  (A:  =  1,2) 

/f  =  0  /i^  =  0 


(9.42) 


and  ak  {k  =  1,2)  are  unknown  coefficients.  Their  values  are  determined 
by  requiring  that  the  corresponding  solution  satisfies  the  desired  boundary 
conditions,  namely  at  i  =  1  and  i  =  AT.  This  requirement 

leads  to  the  following  linear  equation  system: 


—  f^) 


A  O  A 

/\z  2  3Az  2 


(9.43) 


+  "2(1  -  ^/^_3/2  +  ^/^-l/2) 


(9.44) 


which  is  inverted  for  all  Fourier  modes  with  non- vanishing  moduli.  For  the 
zeroth  mode,  i.e.  when  kx  =  ky  =  0,  the  equations  are  integrated  with  the 
boundary  conditions 


dL02.  _  Idp 

dz  ndx 

Sir' = = 0 


(9.45) 

(9.46) 

(9.47) 


A  similar  technique  is  used  to  update  U2  and  oq.  The  discrete  system  of 
equations  is  first  expressed  as: 


=  RzWiQTi) 

i3i5”A'=<3j(c;"+‘) 

^2 


(9.48) 

(9.49) 


with  boundary  conditions 


=  0 


=  0 
if  1^1  ^  0 


=  -^h 

^iN  = -^\n  if  |fc|  7^  0 
^|i  =  0  if|^|  =  0 

^u  =  0  if  1^1  =  0 


I 


Here,  denotes  the  Helmoholtz  operator, 

1 


2ReAz^^^^-^  ^At 


+ 


1 


^  r->n+l 

2ReAz'^ 


ReAz^ 

and  RzW  is  the  corresponding  source  term, 

wc:?)  = 

_ L|m2;j;"+  *  wiJLi  -  2i3;;‘ + g; j;!",' 

2Re''‘'  ““  +  2Jie  Ai? - 

Meanwhile,  is  a  discrete  “Laplacian”  operator  defined  by; 


U2^^l 

lzu2::i  =  *-= 


^2,_.  3 


and 


<53(cni^')  = 


rT^Ti+l  r^n+l 
^It+l  -<^u-i 

Az 


(9.50) 

(9.51) 

(9.52) 


(9.53) 


(9.54) 

(9.55) 

(9.56) 


is  the  corresponding  source  term. 


For  Fourier  modes  having  a  non-vanishing  modulus,  boundary  Green’s 
functions  are  used  to  decouple  the  governing  equations  for  vorticity  and  ve¬ 
locity.  We  first  solve  the  auxiliary  system: 


r  =  R^Wiul’l) 

\  =  011+^  =  0 

1  =  0 

and  then  express  the  solution  of  the  whole  system  as: 

r  +  eLi 

[u2^:i=n,^:i+EU/3,xi, 

Here,  k*  and  are  the  solutions  of  the  elementary  problems: 

'H3K^  =  0  (A:  =  1,2) 

^  /cj  =  1  =  0 

[kI  =  o  4  =  1 

rL3Af  =  0  (A:  =  1,2) 

\  Af  =  0  A^  =  0 


(9.57) 

(9.58) 

(9.59) 

(9.60) 

(9.61) 


and  Pk  {k  =  1,2)  are  unknown  coefficients  whose  values  are  determined  by 
requiring  that  the  corresponding  solution  satisfied  the  desired  wall  boundary 
conditions.  This  requirement  leads  to  the  following  linear  equation  system: 


/?i(l  +  ~  x'^As)  -I- 

Az  2  3Az  2 


3 


^2(1 — As  — 


1 

3A2: 


_ 7/;:”+i  _i_  _JL_j7'”+i 

.  “23  +  o  A  ^2  5 

Az  2  3A2;  2 


+ 


^2(1  —  ■^4-1/2  + 


^  \2 

3a/^ 


(9.62) 

-3/2) 


which  is  inverted  for  every  mode. 


For  the  zeroth  mode,  the  equations  are  subjected  to  the  following  bound¬ 
ary  conditions, 
d(^i. 

=  0  (9.64) 

=  0  (9.65) 

and  inverted  independently  of  the  others. 


The  velocity  compoment  uz  is  updated  by  the  discrete  equation: 


=  iky^u^^  - 

(9.66) 

with  the  boundary  condition 

=  0 

(9.67) 

Finally,  the  temperature  field  is  advanced  by  inverting  the  discrete  elliptic 

equation: 

Gfn+l  ^ 

(9.68) 

subject  to  the  boundary  condition. 

I  f  1,  =  W„(f.  -  1)  if|fc|  =  0 

(9.69) 

II  II 

1  1 

1 

°  II 

o 

(9.70) 
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Here,  G  denotes  the  one-dimensional  Helmholtz  operator, 

_j _ ^ _ '\yn+l  _  _ ^ _ rpn+1 

2RePrAz‘^^  *  2RePrAz‘^ 

RT  the  corresponding  source  term. 


(9.71) 


ufn  _  iL  _  t”  _ Tl"  _  fci+i  J-a-i  i  \K\2ffm 


1  T”  —  n'f'n  I  rpn 

H - ^  H - -t-  $”)  fg  72') 

2RePr  Az^  2Re^  ^  ^  ’ 


9.5  Initial  Conditions 


Simulations  are  performed  in  channel  of  depth  of  0.2mm.  The  periodic 
lengths  in  thestreamwise  and  spanwise  directions  are  0.47rmm  and  0.27rmm 
respectively.  The  initial  velocity  field  is  obtained  by  superposing  on  the 
one-dimensional  streamwise  mean  velocity  profile  (Fig.  9.1)  a  randomly  gen¬ 
erated  homogeneous  3-dimensional  velocity  field.  The  undisturbed  mean 
velocity  profile  is  calculated  using  the  two-layer  model  discussed  in  chapter 
7,  with  the  centerline  velocity  Uc  =  150m/s.  The  purturbation  velocity  field 
has  a  root  mean-square  value  equal  to  0.02Uc.  For  a  kinematic  viscosity 
=  5.0  X  10®m^/s,  the  Reynolds  number  based  on  centerline  velocity  is 
about  6, 000.  In  order  to  avoid  excessive  resolution  requirements,  a  moder¬ 
ate  value  of  the  Prandtl  number,  Pr  =  7,  is  selected.  The  temperature  field 
is  initialized  using  the  mean  temperature  distribution  solution  obtained  in 
Chapter  7  (Fig.  9.2). 


Another  parameter  we  need  to  be  concerned  is  the  value  of  mean  pressure 
gradient,  i.e.  ||,  which  controls  the  intensity  of  forcing  within  the  channel, 
and  therefore  the  mass  flow  rate  within  the  channel.  When  the  undisturbed 
centerline  velocity  is  given,  the  corresponding  friction  velocity  can  be  ob¬ 
tained  numerically  using  the  two-layer  turbulent  velocity  model  within  the 
channel  (Chapter  7).  Then,  the  mean  pressure  gradient  is  determined 
through. 


dp  _  p  2 
dx 

The  mean  pressure  gradient  with  respect  to  y,  is  set  to  zero. 


(9.73) 


9.6  Results  and  Discussion 

The  computations  are  performed  on  a  grid  with  iVi  x  ATj,  x  =  80  x  80  x  160 
mesh  points.  We  have  let  code  run  for  10, 500  time  steps,  which  is  long 
enough  for  the  perturbation  field  to  undergo  significant  deformation.  The 
time  step  At  is  chosen  to  be  O.IA2:  in  normalized  units,  i.e.  the  Courant 
number  is  approximately  0.1.  The  corresponding  physical  time  step  is  about 
0.818ns.  As  mentioned  in  the  previous  section,  wall  heat  transfer  is  modeled 
using  a  heat  transfer  coefficient.  In  the  simulation,  we  have  selected  a  value 
h  =  lOOW/m^  K,  which  is  the  representative  of  forced  air  cooling  at  low 
speed. 


The  evolution  of  the  mean  flow  is  depicted  in  Figs  9.3-9.5,  which  show 
profiles  of  the  mean  streamwise  velocity  mean  spanwise  vorticity  uJ^,  and 
mean  temperature  respectively.  The  profiles  are  obtained  by  averaging  the 
predictions  along  horizontal  planes,  i.e.  planes  parallel  to  the  x-y  plane.  The 
figures  indicate  that  the  velocity,  vorticity  and  thermal  boundary  layers  ad¬ 
just  a  little  during  the  initial  stages  of  the  computations.  Such  an  adjustment 
is  not  surprising  since  the  mean  flow  is  initialized  using  the  approximate  pro¬ 
files  of  a  two-layer  model. 

Spatial  fluctations  of  the  velocity  and  vorticity  fields  are  shown  Figs.  9.6- 
9.11  which  respectively  shown  profiles  of  r.m.s.  values  of  the  fluctating  part  of 
the  streamwise,  spanwise,  and  cross-stream  velocity  components,  and  of  the 
streamwise,  spanwise,  and  cross-stream  vorticity  components.  The  fluctating 
components  are  found  by  subtracting  from  the  raw  data  the  corresponding 
mean  value  on  horizontal  planes,  and  the  averaging  on  the  same  planes.  The 
figures  indicate  that  the  velocity  and  vorticity  fluctations  also  undergo  an 
adjustment  as  time  evolves,  with  an  appreciable  drop  near  the  centerline  of 
the  channel.  This  adjustment  is  to  be  expected  since  the  initial  field  used  in 
the  computations  was  perturbed  using  a  homogeneous  isotropic  field. 

Spatial  fluctuations  of  the  temperature  field  are  first  illustrated  using  tem¬ 
perature  distribution  plotted  along  streamwise  planes,  i.e.  whose  normal  is 
parallel  to  the  streamwise  x-axis.  Figures  9.12-9.14  show  temperatures  con- 


tours  in  the  plane  x  =  27rd  at  times  t  =  2.05/fs,  t  =  5.32/is,  and  t  =  8.59//S, 
respectively.  The  figures  show  that  the  temperature  boundary  layers  at  both 
the  lower  and  upper  walls  spread  slightly  during  the  initial  stages,  and  that 
the  temperature  distribution  is  significantly  more  deformed  on  the  upper 
wall  than  the  lower  one.  The  shape  of  the  deformation  is  indicative  of  the 
presence  of  concentrated  eddies  within  the  boundary  layers.  Comparison  of 
the  temperature  contours  with  velocity  and  vorticity  distribution  generated 
in  a  similar  fashion  (not  shown)  indicated  that  the  non-uniformities  observed 
near  the  top  of  Figs.  9.12-9.14  are  due  to  the  presence  of  streamwise  vortices 
which  are  periodically  generated  within  the  boundary  layer  and  later  ejected 
into  the  ’’free  stream”. 

The  magnitude  of  the  deformation  indicated  by  the  contour  plots  of  the 
temperature  field  leads  us  to  expect  significant  temperature  fluctations  along 
horizontal  planes.  To  quantify  these  deformations,  Fig.  9.15  compares  cross¬ 
stream  profiles  of  the  peak  temperature  and  the  mean  temperature.  The 
results  show  that  the  largest  departures  occur  within  the  boundary  layer 
structure.  In  particular,  for  the  conditions  of  the  simulation,  the  profiles  for 
the  mean  temperature  and  peak  temperature  achieve  their  maxima  at  ap¬ 
proximately  the  same  location  near  the  upper  boundary,  with  temperature 
rises  above  the  ambient  of  lO^K  and  IS^K,  respectively.  Thus,  the  tem¬ 
perature  fluctuations  can  achieve  amplitudes  as  large  as  30%  of  the  mean 
temperature  rise,  and  thus  represent  a  significant  fraction  of  the  total  heat- 


ing  locally  experienced  by  the  fluid. 


The  present  results  should  be  contrated  with  those  of  the  previous  chap¬ 
ter  which,  for  the  case  of  isotropic  turbulence,  reflected  spatial  temperature 
fluctuations  having  very  small  amplitudes.  We  believe  that  this  contrast  is 
due  to  differences  in  the  flow  mechanisms  which  are  associated  with  the  ob¬ 
served  heating.  In  the  isotropic  case,  rapid  mixing  of  the  mixture  effectively 
reduces  the  likelihood  of  occurrence  of  sharp  temperature  gradients.  Such 
mechanisms  are  clearly  lacking  near  solid  boundaries,  where  concentrated 
eddies  may  be  sustained  for  relatively  large  time  periods  in  regions  of  low 
velocity. 

At  any  rate,  the  result  of  the  present  computations  indicate  that  the 
unsteady  near- wall  flow  dynamics  in  transitional/turbulent  flow  regime  can 
have  a  significant  impact  on  shear  heating  phenomena.  In  particular,  esti¬ 
mates  of  peak  temperature  based  on  mean  temperature  profiles  only  should 
generally  be  expected  to  significantly  underpredict  the  risk  of  mixture  igni¬ 
tion. 

At  present,  it  is  unfortunately  not  feasible  to  rely  on  direct  numerical 
simulations  to  extend  the  predictions  to  the  high-Reynolds-number  high- 
Prandtl-number  regime.  It  would  be  possible,  however,  to  obtain  some  con¬ 
servative  estimates  by  extrapolating  low-Reynolds-number  with  the  help  of 


some  reasonable  assumptions,  such  as  Reynolds-number  independence  and 
linear  dependence  on  the  Prandtl  number.  Nonetheless,  it  would  be  crucial 
that  such  extrapolations  be  tested  against  experimental  results,  and  that  ex- 
periements  are  conducted  in  order  to  carefully  examine,  in  particular,  the 
impact  of  near-wall  dynamics  on  shear  heating  of  the  mixture. 
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Figure  9.3:  The  evolution  of  the  mean  streamwise  velocity 
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Figure  9.4:  The  evolution  of  the  mean  spanwise  vorticity 
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Figure  9.5:  The  mean  temperature  evolution 
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Figure  9.7:  The  evolution  of  root  mean-square  spanwise  velocity. 
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Figure  9.8:  The  evolution  of  root  mean-square  cross-stream  velocity. 


Figure  9.10:  The  evolution  of  root-mean-square  spanwise  vorticity. 
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Figure  9.13:  The  contour  plot  of  temperature  distribution  in  y-z  plane(x  = 
27r(i)  att  =  5.32/xs. 


Appendix  A 


This  appendix  discusses  the  numerical  simulation  of  the  parabolized  equa¬ 
tions  of  motion  for  a  constant  property  mixture.  Section  A1  summarizes 
teh  governing  equations  in  cylindrical  {r,z)  coordinates,  while  Section  A. 2 
provides  the  same  system  after  transformation  (4.23)  is  applied.  Numer¬ 
ical  simulation  of  the  unsteady  (steady)  equations  is  discussed  in  Section 
A.3(A.4). 


A.l  Parabolized  Approximation  for  a  Con¬ 
stant  Property  Mixture 


As  mentioned  in  Section  4.3,  the  parabolized  approximation  of  the  equations 
of  motion  is  based  on  ignoring  the  appropriate  streamwise  gradients.  Imple¬ 
mentation  of  the  approximation  reduces  the  system  of  governing  equations 
to: 


duj 

dt 


d  ,  .  d  ,  .  1 


d  ,a;,. 


(A.1) 
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dT  dT  dT 


a, I*. 

R,P,^dr^  rdr^  IL 


(A.2) 

(A.3) 


where 


Vr  = 

V  = 

Z  j.  Qj. 

^  ~  (-3r) 


(A.4) 

(A.5) 


A.2  Stretched  Coordinates 


When  the  cylindrical  coordinate  system  is  stretched  using  the  transformation 
Eq.  (4.23),  the  governing  system  of  equations  given  above  becomes: 

I  C  ^  ^  I  ^  \  ^  ^  du  ^rdu)  w. 

r  8^  ~ 

dT  dT  dT  _  1  ^  dT  CrdT,  , 


where, 


1 

^  r  dz 
ir  ^ 


(A.9) 

(A.10) 


where,  and  ^rr  respectively  denote  d^/dr  and  d'^^/dr^. 


A. 3  Unsteady  Simulation 


As  in  Section  4.2,  a  finite-difference  approach  to  the  simulation  of  the  above 
system  of  equations  is  adopted.  Here,  the  non-linear  convection  terms  are 
treated  using  a  3rd-order  Adams-bashforth  scheme,  while  a  2nd-order  Crank- 
Nicolson  scheme  is  applied  to  the  remaining  terms. 


To  simplify  the  presentation,  we  shall  employ  the  following  abbreviations: 

f 

^  9i  =  (A.  11) 

.  n  =  r(6) 

in  conjunction  with  standard  finite-difference  notation.  With  these  defini¬ 
tions,  time-stepping  is  summarized  as  follows: 

1.  The  vorticity  in  the  domain  interior  is  updated  by  inverting,  at  every 
streamwise  location,  the  discrete  elliptic  equations: 


(A.12) 


with  Dirichlet  boundary  conditions: 


( 


=  0 

nr!, -mil., • 


Here,  H  denotes  the  one-dimensional  Helmholtz  operator, 

fi  f! 


Hu 


'ij  =  [—(— 

^2i2e^2A^ 


+[ 


LJL 

Re2A^ 

n 


+ 


+ 
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2riAe 

1 

2i?er? 
9i 


-I- 


Ae2 


At 


-F 


fi 


_  _ 

^2Re^A^^  ^  2Ae  ^  2riA^ 


(A.13) 


(A.14) 
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RW  the  corresponding  source  term, 


1 


=  -^  +  :^(235y"-16S»ti  +  55?] 


At  12' - «  '■> 

£  ,  •T1-+2  .  ,71+2 

+te  +  r)-^±M7XT^-^ 


(A.15) 


_  (wiie)i+ij  -  (wtifjj-ij  ,  (w»^)i+ij  -  j 

2Ae  +  2Ai - 

the  vorticity  convection  term. 


(A.16) 


2.  The  new  streamfunction  distribution  is  computed  by  inverting,  at  every 
streamwise  section,  the  discrete  “Laplacian”  opeartors: 


=  r,a;rt3 


(A.17) 


with  the  Dirichlet  boundary  conditions  given  in  Eqs.  (4.6)  amd  (4.7). 
Here, 


Ub  =  -  Si!)  _  u-  2//  / 

I  /  /»  2/f 

4riA^  2A^ 

is  the  discrete  Laplacian  associated  with  Eq.  (A-7). 

3.  Determine  the  solid  surface  vorticity  values  using: 

^.n+3  _  o  J-2  ~  i’N'l-l, 

‘^Nr,i  -  -^JNr - T7^ - 


(A.18) 


(A.19) 


4.  Update  the  velocity  field  using: 


{ 


2riAz 

(Vr\ti  =  fi'ti±LiZjb=ld. 
\^z)tj  Jt  2rjA5 


(A.20) 


5.  Advance  the  temperature  distribution  by  inverting,  at  every  streamwise 
location,  the  discrete  elliptic  equations: 


(A.21) 


subject  to  the  Neumann  conditions  given  in  Eqs.  (4.6-4.7).  Here,  G 
denotes  the  one-dimensional  Helmholtz  opearator. 
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'^^RePr  A^^ 
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RT  the  corresponding  source  term, 
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2RePr^-'^  Ae2 

r  rpn-\-2  rpn’j-2  771 


(A.22) 


+  (Pi  +  ffl— (A.23) 


and, 


Pli  = 


2Az 


(A.24) 


the  thermal  convection  term.  Note  that  the  viscous  dissipation  function 
is  approximated  through: 


I  ^2  ^  (A.25) 

while  the  grid  stretching  derivatives  /j  and  Qi  are  evaluated  numerically, 
respectively  using: 


/,= 

’"t+l  ~  ^t-1 

f, _ /^i+1  ~  2ri  -h  ri_i .  /ffi+l  ~  I'i-l  >^3 

)/(  2A5  > 

(A.26) 

(A.27) 

A.*4  St©£icly  Cod© 

When  considering  the  steady  parabolized  equations  of  motion,  the  streamwise 
coordinate  2:  may  be  regarded  as  a  time  variable  and  the  solution  is  marched 
fro  mthe  inlet  of  nozzle  towards  the  exit  in  small  Az  increments.  Simulation 
of  the  steady  parabolized  equations,  which  also  employs  the  stretched  grid 
technique  discussed  above,  on  the  following  finite-difference  discretizations 
of  the  vorticity  transport,  streamfunction  and  energy  equations: 

(A.28) 

t  M77‘+'  =  +  f  A#?+‘ 

where 
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K)»-+i  Lz.  fl 


2  Ae  i2e  Ae2^2Ae^2riAe^^ 


^2A^  A^^  2riAe^^*“^^  A^”2^‘ 

'^^2riAe  “  ^  ~  A^]^»+i 


(A.30) 


jj^jin+l  ^  r  r  Az  /j2  gj  /.  ^ 

‘  ^  2  Ae  RePr^Ae  2Ae  2riAe^J 

2  Ae  ii:ePr^Ae2^2Ae  2riArJ  '+1 
and  the  notation  w”  =  uj{^i,Zn)  =  w(iA^,nAz)  is  used.  Meanwhile,  the 
velocity  components  v^  and  v^  appearing  in  Eqs.  (A29-31)  are  respectively 
approximated  using: 


n  oA^+1 


-‘fp 


"  ViAz 

(,. «+.  _ ,  c-S‘  -  e*i‘ 
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(A.32) 

(A.33) 


As  in  the  unsteady  case,  solution  marching  requires  the  inversion  of  three 
elliptic  operators  for  the  vorticity,  streamfunction  and  temperature.  How¬ 
ever,  in  the  present  steady  case,  implicit  vorticity  boundary  conditions  are 
implemented,  so  that  the  equation  systems  for  vorticity  and  streamfunction 
are  coupled  at  the  boundary.  In  order  to  reduce  the  computational  cost  as¬ 
sociated  with  simultaneous  inversion  of  the  two  systems,  a  boundary  Green’s 
function  technique  is  incorporated  into  the  computations.  This  technique 
exploits  the  linearity  of  the  elliptic  operators  by  first  solving  the  auxiliary 
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systems: 


\  = uif = 0 

'  = co^+\ 

<  =  0 

^  =foru^(r)dr 

and  then  expressing  the  solution  to  the  full  system  as: 
r  a;r+i  =  a;”+'  +  6/Ci 

1  +  eAi 

where  /Cj  and  Aj  are  the  solutions  to  the  elementary  problems: 

Hki  =  0 
<  Ki  =0 
.  «Arr  =  1 

(  LXi  =  KiTi 

\  «1  =  «Wr  =  0 


(A.34) 

(A.35) 

(A.36) 

(A.37) 

(A.38) 


and  e  is  an  arbitrary  constant.  Choosing  e  so  that  vorticity  boundary  condi¬ 
tions  are  satisfied,  we  find: 


_  ,7,n+l 

_  Agii  \yNr  ^Nr-X) 


e  = 


1  A 

■I  +  -^'^Nr-1 


(A.39) 


No  such  difficulty  is  encountered  in  the  marching  of  the  energy  equation,  since 
homogeneous  Neumann  boundary  conditions  are  imposed  both  at  r  =  0  and 


r  =  1. 


Appendix  B 


This  appendix  discusses  the  numerical  simulation  if  the  parabolized  equations 
of  motion  for  an  incompressible  variable  property  mixture.  Two  models  are 
used  in  the  computations:  the  first  model  accommodates  a  temperature- 
dependent  viscosity,  but  assume  that  a  constant  thermal  conductivity.  In 
the  second  model,  both  the  viscosity  and  thermal  conductivity  are  taken  to 
be  temperature  dependent. 

Since  the  former  model  may  be  regarded  as  a  special  case  of  the  second,  only 
the  more  general  case  is  discussed.  Section  B1  summarizes  the  governing 
equations  for  the  parabolized  approximation  in  cylindrical  (r,  z)  coordinates, 
while  Section  B2  provides  the  same  system  after  transformation  (4.23)  is  ap¬ 
plied.  Numerical  simulation  of  the  unsteady  (steady)  equations  is  discussed 
in  Section  B3  (B4). 


B.l 


Parabolized  Approximation  for  a  Vari¬ 
able  Property  Mixture 


When  the  appropriate  streamwise  gradients  are  ignored,  the  governing  equa¬ 
tions  given  in  Section  5  reduce  to: 
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(B.l) 

(B.2) 

(B.3) 

(B.4) 

(B.5) 

(B.6) 

(B.7) 


and  the  is  used  to  denote  a  dimensional  quantity. 


B.2  Stretched  Coordinates 


When  the  cylindrical  coordinate  system  is  stretched  using  the  transformation 
Eq.  (4.23),  the  governing  system  of  equations  given  above  becomes: 


r  dC  r2^  ^ r’ 


2^^^*  *  du 


.2  52^  .  V.  dip 


+  ^r 
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where 


(B.8) 

(B.9) 

(B.IO) 


»« =  -if 

—  S^Si. 

^  r 


$ 


while  and  f..  respectively  denote  d^/dr  and  di^^/dr^. 


(B.ll) 

(B.12) 


B.3  Unsteady  Simulation 

As  in  Appendix  A,  a  finite-difference  methodology  is  adopted  in  the  simula¬ 
tion  of  the  above  system  of  equations.  Here,  teh  non-linear  convection  terms 
are  treated  using  a  3rd-order  Adams-Bashforth  scheme,  while  a  2nd-order 
Crank-Nicolson  scheme  is  applied  to  the  remaining  terms.  Using  the  same 


notation  of  Appendix  A,  time-stepping  is  summarized  as  follows: 


1.  The  vorticity  in  the  domain  interior  is  updated  by  inverting,  at  every 
streamwise  location,  the  discrete  elliptic  equations: 


with  Dirichlet  boundary  conditions: 

/  =  o 

[  <3  = 

Here,  H  denotes  the  one-dimensional  Helmholtz  operator, 
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RW  the  corresponding  source  term, 
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(B.16) 


and, 


-  Ji  2Ae 
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{u)Vz)iJ+i  {u}Vz)ij^i 

2Az 


(B.17) 


the  vorticity  convection  term. 


Note  that  the  viscosity  values  at  time  level  n  +  3,  which  appear  in  the 
definition  of  the  Helmholtz  operator  (B.15),  are  not  generally  known 
before  the  energy  equation  is  solved.  In  order  to  avoid  a  nonlinear 
coupling  between  the  vorticity  transport  and  energy  equations,  the  fol¬ 
lowing  approximation  is  used  in  the  computations: 


{«•)”+’  =  2(v-)r+^  -  («•)”+' 


(B.18) 


2.  The  new  streamfunction  distribution  is  computed  by  inverting,  at  every 
streamwise  section,  the  discrete  “Laplacian”  operator: 


= 


,71+3 

hJ 


(B.19) 


174 


with  the  Dirichlet  boundary  conditions  given  in  Eqs.  (4.6)  and  (4.7). 


LUJ  ■  =  (Ji-  -  ^  I 

4-  (  9i  Vis, 

^2riA^  2A^ 

is  the  discrete  Laplacian  associated  with  Eq.  (B.9). 

3.  Determine  the  solid  surface  vorticity  values  using: 


n+3  _  _0  f2  ‘^N^d  ~ 
JNr 


4.  Update  the  velocity  field  using: 


—  2riAz 
K<Jz)t,j  —  Jt  2r.Af 


(B.20) 


(B.21) 


(B.22) 


5.  Advance  the  temperature  distribution  by  inverting,  at  every  streamwise 
location,  the  discrete  elliptic  equations: 


GTr^f^  = 


(B.23) 


subject  to  the  Neumann  conditions  given  in  Eqs.  (4.6-4. 7).  Here,  G 
denotes  the  one-dimensional  Helmholtz  operator, 
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and, 


2A^ 

the  thermal  convection  term 
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(B.25) 


(B.26) 


Note  that  the  viscosity  and  thermal  conductivity  values  at  time  level 
n  +  3,  which  appear  in  the  definition  of  the  Helmholtz  opeator  (B.24) 
and  the  cource  term  (B.25),  are  not  known  before  the  new  temperature 
distribution  is  determined.  In  order  to  avoid  nonlinear  iterations,  the 
following  approximations  is  used  in  the  computations: 


(r)^t3  =  2{k*)if  - 


(B.27) 

(B.28) 


As  in  Appendix  A,  the  viscous  dissipation  function  is  approximated 
through: 


while  the  grid  stretching  derivatives  /,•  and  Qi  are  evaluated  numerically, 
respectively  using: 
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(rj+i  -  2ri  +  ri-i)/A^^ 
[(ri+i  -  ri_i)/2A^]3 


(B.30) 

(B.31) 


B.4  Steady  Code 


Simulation  of  the  steady  parabolized  equations  of  motion  follows  a  sim¬ 
ilar  approach  to  that  used  for  constant  property  flow.  Specifically, 
the  simulation  is  based  on  the  following  stretched  grid  finite  difference 
approximations  of  the  vorticity  transport,  streamfunction  and  energy 
equations: 
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and  the  notation  w”  —  w(^i,  2:„)  =  ^(iA^,  nAz)  is  used.  The  velocity 
components  and  appearing  in  Eqs.  (B.33-35)  are  respectively 
calculated  using: 


n+l  _ 


Ti^Z 


n+1 


2riAe 


(B.36) 

(B.37) 


while  the  unknown  values  (1/*)"+^  and  appearing  in  the  defini¬ 

tions  (B.33)  and  (B.35)  are  respectively  approximated  through: 


{u*)r^  =  2{u*Y,  -  {u*rr^  (B.38) 

(A:*)r'  =  2(r)^  -  (r)r^  (B.39) 


The  numerical  solution  of  Eqs.  (B.33-35)  uses  the  same  LU-decomposition 
in  conjunction  with  Boundary  Green’s  function  technique  summarized 
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in  Appendix  A4. 
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